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

K Liu, A Parthasarathy, M Keith, C Tiburzi, S.C Susarla, J Antoniadis, A Chalumeau, S Chen, Ismaël Cognard, A Golden, et al.

► To cite this version:

K Liu, A Parthasarathy, M Keith, C Tiburzi, S.C Susarla, et al. The impact on astrometry by solar-wind effect in pulsar timing. Monthly Notices of the Royal Astronomical Society, 2025, 536 (3), pp.2603-2617. 10.1093/mnras/stae2727. hal-04891923

HAL Id: hal-04891923

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

Submitted on 21 Mar 2025

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

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

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The impact on astrometry by solar-wind effect in pulsar timing

K. Liu<sup>1,2*</sup> A. Parthasarathy<sup>3,2</sup> M. Keith<sup>4</sup> C. Tiburzi<sup>5</sup> S. C. Susarla<sup>6</sup> J. Antoniadis<sup>7,2</sup>

A. Chalumeau<sup>3</sup>, S. Chen<sup>1</sup>, I. Cognard<sup>8</sup>, A. Golden<sup>6</sup>, J.-M. Grießmeier<sup>8,9</sup> L. Guillemot<sup>9</sup>

G. H. Janssen<sup>4,10</sup> E. F. Keane<sup>11</sup> M. Kramer<sup>2,4</sup> J. W. McKee<sup>12</sup> M. B. Mickaliger<sup>4</sup>

G. Theureau<sup>8,9,13</sup> and J. Wang<sup>14</sup>

Accepted 2024 December 6. Received 2024 December 5; in original form 2024 August 19

## ABSTRACT

Astrometry of pulsars, particularly their distances, serves as a critical input for various astrophysical experiments using pulsars. Pulsar timing is a primary approach for determining a pulsar's position, parallax, and distance. In this paper, we explore the influence of the solar wind on astrometric measurements obtained through pulsar timing, focusing on its potential to affect the accuracy of these parameters. Using both theoretical calculation and mock-data simulations, we demonstrate a significant correlation between the pulsar position, annual parallax, and the solar-wind density parameters. This correlation strongly depends on the pulsar's ecliptic latitude. We show that fixing solar-wind density to an arbitrary value in the timing analysis can introduce significant bias in the estimated pulsar position and parallax, and its significance is highly dependent on the ecliptic latitude of the pulsar and the timing precision of the data. For pulsars with favourable ecliptic latitude and timing precision, the astrometric and solar-wind parameters can be measured jointly with other timing parameters using single-frequency data. The parameter correlation can be mitigated by using multifrequency data, which also significantly improves the measurement precision of these parameters; this is particularly important for pulsars at a medium or high ecliptic latitude. Additionally, for a selection of pulsars we reprocess their European Pulsar Timing Array Data Release 2 data to include modelling of solar-wind effect in the timing analysis. This delivers significant measurements of both parallax and solar-wind density, the latter of which are consistent with those obtained at low-frequency band. In the future, combining pulsar timing data at gigahertz and lower frequencies will probably deliver the most robust and precise measurements of astrometry and solar wind properties in pulsar timing.

**Key words:** astrometry – solar wind – pulsars: general.

## 1 INTRODUCTION

Pulsars are highly magnetized, 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; Guo et al. 2021; Kramer 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

<footnote n="1"><sup>1</sup>Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China</footnote>

<footnote n="2"><sup>2</sup>Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote n="3"><sup>3</sup>ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, NL-7991 PD Dwingeloo, the Netherlands</footnote>

<footnote n="4"><sup>4</sup>Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="5"><sup>5</sup>INAF – Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius (CA), Italy</footnote>

<footnote n="6"><sup>6</sup>Physics, School of Natural Sciences, Ollscoil na Gaillimhe – University of Galway, University Road, Galway H91 TK33, Ireland</footnote>

<footnote n="7"><sup>7</sup>FORTH Institute of Astrophysics, N. Plastira 100, 70013 Heraklion, Greece</footnote>

<footnote n="8"><sup>8</sup>Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans/CNRS, F-45071 Orléans Cedex 02, France</footnote>

<footnote n="9"><sup>9</sup>Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, F-18330 Nançay, France</footnote>

<footnote n="10"><sup>10</sup>Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, NL-6500 GL Nijmegen, the Netherlands</footnote>

<footnote n="11"><sup>11</sup>School of Physics, Trinity College Dublin, College Green, Dublin 2, D02 PN40, Ireland</footnote>

<footnote n="12"><sup>12</sup>Physics and Astronomy Department, Union College, 807 Union Street, Schenectady, NY 12308, USA</footnote>

<footnote n="13"><sup>13</sup>Laboratoire Univers et Théories LUTH, Observatoire de Paris, Université PSL, CNRS, Université de Paris, F-92190 Meudon, France</footnote>

<footnote n="14"><sup>14</sup>Ruhr-Universität Bochum, Fakultät für Physik und Astronomie, Astronomisches Institut (AIRUB), D-44801 Bochum, Germany</footnote>

<footnote>*E-mail: kliu.psr@gmail.com</footnote>
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, Manchester & Wang 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, Chatterjee & Cordes 2013).

Typically, there are two primary ways to precisely measure pulsar astrometry. The first method involves the canonical very-long-baseline-interferometer (VLBI) technique that delivers precise localization 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 milliarcsecond (mas; Ding et al. 2023).

The other method involves modelling the annual parallax of the pulsar via precision pulsar timing (e.g. EPTA Collaboration 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, Hobbs & Manchester 2006). For a canonical pulsar, its position can be measured with a precision of 10–100 mas (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 microarcsecond (µas; e.g. EPTA Collaboration 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 data set. 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, Tiburzi & Zucca 2020; Tiburzi et al. 2021). Typically, in pulsar timing the solar-wind effect is described by a spherical model parametrized by the number density of free electrons at a 1 au distance from the Sun ($n_{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_{sw}$ was used as a fixed constant input (e.g. Reardon et al. 2021; EPTA Collaboration 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^{-3}$, while in the TEMPO2 software package it is $4 \times 10^{-3}$. Using the NANOGrav 11-yr data set, Madison et al. (2019) suggested an average value of $7.9 \times 10^{-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_{sw}$ value to different pulsars and data sets is probably not optimal. Moreover, as we point out later, if a fixed constant $n_{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 European Pulsar Timing Array (EPTA) ‘DR2full’ data set guided by our findings is shown in Section 4. We conclude in Section 5 with a brief discussion.

## 2 THE CORRELATION BETWEEN ASTROMETRIC AND SOLAR-WIND PARAMETERS

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

$$T_p = t_{\text{topo}} - \Delta_{\text{R}\odot} - \Delta_p - \Delta_{\text{sol}} - \dots, \quad (1)$$

where $T_p$ is the pulsar proper time and $t_{\text{topo}}$ is the topocentric time. Here $\Delta_{\text{R}\odot}$, $\Delta_p$, and $\Delta_{\text{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_{R\odot}$, can be written as

$$ \Delta_{R\odot} = -\frac{|\vec{r}(E)| \cos \theta(E)}{c}, \qquad (2) $$

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

$$ |\vec{r}(E)| = a(1 - e \cos E), \qquad (3) $$

$$ \cos \theta(E) = \cos \beta \cos(A_T - \lambda), \qquad (4) $$

where $a$, $e$, and $A_T$ are the semimajor 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_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)

$$ \Delta_p = \frac{|\vec{r}(E)|^2 \sin^2 \theta(E)}{2c} \omega. \qquad (5) $$

Here, the parallax $\omega$ 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_{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):

$$ \Delta_{sol} = \frac{DM_{\odot}}{D_0} \cdot \frac{1}{f^2} = \frac{n_{sw}(1 \text{ AU})^2 \theta(E)}{D_0 |\vec{r}(E)| \sin \theta(E)} \cdot \frac{1}{f^2}, \qquad (6) $$

where $D_0 = 2.410 \times 10^{-16} \text{ cm}^{-3} \text{ pc}$ and $f$ is the observing radio frequency. The only model parameter here, $n_{sw}$, is the mean electron density at 1 au heliocentric radius. Though in reality, the solar wind exhibits complex time variability (e.g. Sokól et al. 2013; Provornikova et al. 2014; Shaifullah et al. 2020; Tiburzi et al. 2021; Susarla et al. 2024), 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. Fig. 1 gives the waveform of $\Delta_p$ and $\Delta_{sol}$ over the course of a year, for various ecliptic latitudes of the pulsar. It can be seen that the two waveforms follow anticorrelated 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

$$ \begin{aligned} \rho_{\Delta_p, \Delta_{sol}} &= \frac{\overline{\Delta_p \cdot \Delta_{sol}} - \overline{\Delta_p} \cdot \overline{\Delta_{sol}}}{\sqrt{\overline{\Delta^2_p} - (\overline{\Delta_p})^2} \sqrt{\overline{\Delta^2_{sol}} - (\overline{\Delta_{sol}})^2}} \\ &= \frac{\overline{P \cdot Q} - \overline{P} \cdot \overline{Q}}{\sqrt{\overline{P^2} - (\overline{P})^2} \sqrt{\overline{Q^2} - (\overline{Q})^2}}, \end{aligned} \qquad (7) $$

where

$$ \Delta_p = P \cdot \omega, \quad \Delta_{sol} = Q \cdot n_{sw}, \qquad (8) $$

and

$$ P = \frac{|\vec{r}(E)|^2 \sin^2 \theta(E)}{2c}, \quad Q = \frac{(1 \text{ AU})^2 \theta(E)}{D_0 |\vec{r}(E)| \sin \theta(E)} \cdot \frac{1}{f^2}. \qquad (9) $$

<figure><img src="image_2.png" /><figcaption>Figure 1. Time delay in pulsar timing caused by annual parallax ($\Delta_p$, in solid lines) and solar wind effect ($\Delta_{sol}$, in dashed lines), as a function of days after the perihelion of the Earth. A few values of the pulsar ecliptic latitude $\beta$ were applied, separately, while the ecliptic longitude was set to be $\lambda = 0$ deg.</figcaption></figure>

The overline in equation (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

$$ \overline{F(E)} = \int_{0}^{2\pi} F(E)(1 - e \cos E) dE. \qquad (10) $$

As shown in equation (1), the timing model is a linear combination of a sequence of delay terms including $\Delta_p$, $\Delta_{sol}$. Note that the two model parameters, $\omega$ and $n_{sw}$, are simply the amplitude of the two terms, $\rho_{\Delta_p, \Delta_{sol}}$ should directly reflect the correlation of the two parameters in a least-square fit to the timing model. Meanwhile, since $\Delta_p$ and $\Delta_{sol}$ follow an anticorrelated trend, an increase in $\Delta_p$ (thus in $\omega$) requires a decrease in $\Delta_{sol}$ (thus increase in $n_{sw}$) to keep $\chi^2$ around its minimum in a least-square fit. This means that $\omega$ and $n_{sw}$ are positively correlated, i.e. $\rho_{\omega, n_{sw}} = -\rho_{\Delta_p, \Delta_{sol}}$. Thus, it is clear from equation (7) that $\rho_{\omega, n_{sw}}$ does not depend on the value of $n_{sw}$ and $\omega$ 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_{\omega, n_{sw}}$ is likely to exhibit a significant dependency on $\beta$.

Note that the correlation estimate above only takes into account $\Delta_p$ and $\Delta_{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 $\omega$ and $n_{sw}$, the $\rho_{\omega, n_{sw}}$ calculated from the timing analysis can be significantly different. None the less, 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 equation (A12). The derivation including only $\omega$ and $n_{sw}$ gives equation (A16), exactly the same as in equation (7).

Fig. 2 shows the correlation coefficients $\rho_{\omega, n_{sw}}$, $\rho_{\lambda, n_{sw}}$, $\rho_{\beta, n_{sw}}$ calculated from equations (A12)–(A14) and equation (A16) for different ecliptic latitudes of the pulsar. It can be seen that when
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Correlation coefficient of astrometric and solar-wind model parameters, all as a function of ecliptic latitude of the pulsar. The calculations were carried out based on equations (A12) and (A16). Here, the ecliptic longitude $\lambda$ was set to 45 deg.</figcaption></figure>

the position parameters are also included in the fit, $\rho_{\sigma\tau, n_{sw}}$ is indeed significantly altered; it becomes larger for all $\beta$ values shown in the figure. The correlation between $\beta$ and $n_{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. First, the estimates above assume single-frequency observations, whereas in reality, multifrequency 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.

## 3 SIMULATION ANALYSIS

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 data set also encompasses a set of other standard timing parameters that 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 multifrequency 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 that 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.

### 3.1 Correlation on more realistic occasions

In a global timing fit alongside other model parameters, the correlations between position, parallax and solar-wind parameters are expected to align closely with the derivations in equations (A12)–(A14), as correlations with other parameters are generally minial. Here, we evaluated the correlations between astrometric and solar-wind parameters under more realistic timing-analysis scenarios. The standard timing model as described above was used to simulate the data at L-band frequency (1.4 GHz), with high cadence (ten TOAs per day, randomly distributed in time) and precision (1 ns) to fully capture the sharp feature of solar-wind delay during solar conjunction and accurately estimate the correlation coefficients. The results are presented in Fig. 3. It can be seen that $\rho_{\sigma\tau, n_{sw}}$ obtained from the simulation is almost identical as that from equation (A12). Regarding the position parameters, $\rho_{\lambda, n_{sw}}$, $\rho_{\beta, n_{sw}}$ also exhibit a very similar trend as in Fig. 2. Meanwhile, it should be noted that here the values of $\rho_{\beta, n_{sw}}$ are generally lower than those from equation (A14) at middle and low ecliptic latitudes ($\beta \lesssim 50$ deg). This is mainly because the positional variation in time (i.e. proper motion) mitigates the degeneracy between the position and solar-wind parameters. The proper motion parameters, $\mu_{\lambda}$, $\mu_{\beta}$ start to show significant correlation with $n_{sw}$ only when $\beta$ is approaching 70 deg. In the case where the equatorial coordinate system is used in the timing model, both the position parameters, namely right ascension ($\alpha$) and declination ($\delta$), become significantly correlated with $n_{sw}$ when $\beta$ is larger than a few degrees. Meanwhile, $\rho_{\sigma\tau, n_{sw}}$ follows exactly the same trend as in the ecliptic coordinate system, as expected to be independent of the chosen coordinate system. For the remainder of this paper, we will adhere to the ecliptic coordinate system.

As anticipated from equations (A12)–(A14), the correlations in the case of a single observing frequency as above are independent of choice of the frequency and timing precision of the data. Meanwhile, since the solar-wind delay is chromatic while astrometric terms are not, utilizing data with multifrequency coverage can effectively mitigate the correlation between astrometric and solar-wind parameters. The extent of this mitigation depends on the effective frequency coverage of the data which is determined by the relative frequency separation and precision of the data at different frequency bands. For a detailed investigation, we conducted simulations considering a dual-band scenario. In each iteration, the central frequency and timing precision of the reference band were fixed to 1.4 GHz and 200 ns, respectively, while those of the other band were set to a different pair of values. Then the corresponding $\rho_{\sigma\tau, n_{sw}}$ were obtained by the timing fit using this dual-band data set. Note that as indicated by equation (A17), $\rho_{\sigma\tau, n_{sw}}$ is decided by the ratio of the frequency and timing precision at these two bands but not their
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Correlation coefficient of $n_{sw}$ and other astrometric parameters, as a function of ecliptic latitude of the pulsar. The top plot used the ecliptic coordinate system in the timing model, while the bottom employed equatorial coordinate system. The ecliptic longitude was set to 45 deg in these simulations.</figcaption></figure>

absolute value, the results of these simulations are also applicable to other choices of the frequency and timing precision of the reference band.

Fig. 4 shows the $\rho_{\phi,n_{sw}}$ for different ratios of the frequency and precision of the alternative band to those of the reference band. The results indicate that for very small $\beta$, the correlation is generally low ($< 0.25$) in all cases. However, as the ecliptic latitude increases, effective multifrequency coverage becomes much more beneficial in mitigating the correlation. For instance, at $\beta = 15$ deg, it can decrease the correlation coefficient from above 0.7 to below 0.3. At 45 deg, the reduction is even substantial, from over 0.9 down to 0.1. Additionally, for a given frequency ratio, the correlation coefficient is lowest when the high-frequency band has a relatively higher

timing precision. This finding suggests that achieving better timing precision at higher frequencies is crucial for effectively mitigating the correlation between the astrometric and solar-wind parameters in multifrequency pulsar timing analysis. It also aligns with the results from Lee et al. (2014); Donner et al. (2020) showing the importance of obtaining higher timing precision in high-frequency observations for DM correction, which requires more dedicated efforts compared to those at low frequencies.

## 3.2 Bias in astrometry for fixed $n_{sw}$

In timing analysis, it is common to, when applying the spherical model as the correction for the solar-wind delay, give a fixed free-electron density value for the solar wind (e.g. Edwards et al. 2006; Reardon et al. 2021; EPTA Collaboration 2023a). If this assumption differs from the actual time-averaged value, it would introduce bias in other fitted timing parameters, primarily the astrometric parameters, due to their significant correlations. The bias could be substantial compared to the measurement uncertainty of the parameters, depending on the amplitude of the solar-wind effect primarily decided by the observing frequency and ecliptic latitude of the pulsar compared to the timing precision of the data. To investigate this issue, we conducted a series of simulations with various choices of ecliptic latitude of the pulsar, offset in $n_{sw}$ from the true value, frequency and precision of the TOAs. The solar-wind density was assumed to be time-invariant across the data set.

Fig. 5 illustrates the absorption of the solar-wind delay feature in the timing data by astrometric parameters. It shows the pre-fit and post-fit timing residuals where $n_{sw}$ was fixed to $6 \text{ cm}^{-3}$, $2 \text{ cm}^{-3}$ off from the value ($8 \text{ cm}^{-3}$) used to simulate the data. It can be seen that for a low ecliptic latitude of $\beta = 3$ deg, the feature is mostly retained after the fit as the correlation is moderate as shown in Fig. 3. When $\beta$ increases, the correlation grows rapidly, and thus the remaining feature shrinks after absorption by other timing parameters. For $\beta = 15$ deg, the post-fit residuals are around the 100-ns scale, about half of the pre-fit residuals scale. Note that 100 ns is approximately the best timing precision achieved so far (e.g. Liu et al. 2020), suggesting that in this case even with single-frequency data the solar-wind feature could still be measurable with the current best timing precision. Meanwhile, for high ecliptic latitude, e.g. $\beta = 45$ deg, only a small fraction of the solar-wind feature remains after the fit for other timing parameters. This indicates that it would be difficult to break the degeneracy between parallax and solar-wind parameters with single-frequency data in this case.

Fig. 6 shows the bias in the astrometry parameters caused by the attempt to model the feature caused by various $n_{sw}$ offset values ($\Delta n_{sw}$) with astrometric parameters. It can be seen that as the ecliptic latitude increases, reducing the solar-wind delay, the bias decreases. The bias in ecliptic longitude is generally on the order of 1 µas, smaller than that in ecliptic latitude and parallax, as expected from its lower correlation with $n_{sw}$. The bias in ecliptic latitude is a few orders of magnitude larger, while strongly depending on the value of ecliptic latitude itself. The bias in parallax is generally a few times of 0.1 mas, which can be critical for some precisely timed pulsars whose timing parallax measurement can be more precise than 0.1 mas.

The bias can be either significant or negligible compared to the measurement error of the parallax, strongly depending on the ecliptic latitude of the pulsar, the precision of the data set, and the deviation of $n_{sw}$. Fig. 7 shows the normalized bias in parallax by its measurement error for various timing precisions and offsets in $n_{sw}$, across a few values of $\beta$. Here, for each simulation run, we generated a 5-yr
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Correlation coefficient of $n_{sw}$ and $\omega$ from data with dual-frequency coverage, for various frequency and residual rms ratios of the two bands. These of the reference band were set to 1.4 GHz and 200 ns, respectively, in the simulations. The secondary Y-axis shows the frequencies of the alternative band. However, since the correlation coefficient only depends on the frequency ratio of the two bands, these values are directly scalable for other options of the reference frequency. Note that the colour scale in all the panels are identical, which deepens as the value increases.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Pre-fit (blue circles) and post-fit (black dots) timing residuals, showing the timing signature of the solar wind from simulations with different ecliptic latitude of the pulsar. Here, the TOAs were simulated at 1.4 GHz and with $n_{sw} = 8 \text{ cm}^{-3}$, while the residuals were created by fixing $n_{sw}$ to $6 \text{ cm}^{-3}$ but fitting all the other timing parameters. The ecliptic longitude here was set to 45 deg. The TOA precision was assumed to be 1 ns in this illustration.</figcaption></figure>

observation with two TOAs per week at 1.4 GHz. It can be seen that when the ecliptic latitude is as low as a few degrees, the bias can be significant, e.g. above $3\sigma$, for a timing precision around 1 µs or even above. At $\beta = 15$, it is significant only when the timing precision is a few hundred nanoseconds or better. For $\beta$ at 45 deg, the bias is generally small compared to the measurement error, unless the timing precision is better than a couple of hundred nanoseconds

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Bias in ecliptic longitude ($\Delta\lambda$), ecliptic latitude ($\Delta\beta$), and parallax ($\Delta\omega$), induced by the difference in solar-wind density applied in the timing analysis from the actual value ($\Delta n_{sw}$). Here, the results are shown for pulsars of different ecliptic latitudes, and the parallax was set to be 1 mas. The observing frequency of the simulated data was 1.4 GHz. To accurately estimate these bias values, we applied a very high TOA precision (1 ns) in the simulation, although these results should be independent of the assumed timing precision of the data as far as the astrometric parameters are significantly measured. The same applies to the results in Fig. 8.</figcaption></figure>

and the offset in $n_{sw}$ is beyond a few units. In all these cases, the bias significance in ecliptic latitude closely matches that in parallax, as expected from their very similar correlation coefficient shown in Fig. 3. Meanwhile, the bias in ecliptic longitude is always well below $1\sigma$, owing to its weak correlation with $n_{sw}$.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

**Figure 7.** Heat map showing the absolute value of the bias in fitted parallax normalized by its measurement uncertainty for a few $\beta$ values. The colour scales are in all plots are identical so they can be contrasted. Each simulation run here included a 5-yr-long observation with two TOAs per week. The observing frequency of the simulated data set was 1.4 GHz. Note that the colour scale in all the panels are identical, which deepens as the value increases.

The scale of the solar wind delay depends on the observing frequency, suggesting that the bias is frequency-dependent as well. For demonstration, in Fig. 8 we simulated data set within a wide range of observing frequencies (700 MHz to 4 GHz) to inspect how the bias of the astrometric parameter changes. Clearly, for a given offset in $n_{sw}$, the bias in all three astrometric parameters reduces as the frequency increases and the solar-wind effect becomes smaller. Similar to Fig. 7, the results again show that for a small ecliptic

**Figure 8.** Bias in astrometric parameters for given offset in $n_{sw}$ and observing frequency of the data set. The biases in $\lambda$ is in unit of $\mu$s, while the biases in $\beta$ and $\omega$ are in unit of mas. For these simulations, the assumed $\beta$ value is 3 deg and $\lambda$ is 45 deg.

latitude (e.g. $\beta$ = 3 deg) the bias in parallax can be considerable even for observing frequency above 1 GHz. For observations containing data at multiple frequencies, the overall bias typically falls between the biases at the individual frequencies. The exact value depends on the relative timing precision of each data set.

Note that bias analysis above may represent a worse-case scenario. In reality, including modelling of DM variation in the timing analysis could be able to absorb some fraction of the solar-wind feature, thereby reducing the bias in astrometric parameters when a fixed $n_{sw}$ is used. This though would highly depend on the assumed model
<figure><img src="image_14.png" /><figcaption>**Figure 9.** Fractional measurement error of $\bar{\omega}$ and $n_{sw}$ as a function of ecliptic latitude $\beta$. The black lines represent the errors of $\bar{\omega}$ and the blue lines stand for those of $n_{sw}$. The data were simulated at a frequency of 1.4 GHz.</figcaption></figure>

for the DM variation. For instance, the piece-wise DM modelling, commonly referred to as ‘DMX’, is expected to capture DM variation induced by effects including the solar wind, down to a time-scale of a few weeks (e.g. Keith et al. 2013; Agazie et al. 2023). Meanwhile, it is not anticipated that a power-law DM noise modelling will fully absorb this feature in the data, since it exhibits an annual periodicity rather than being a stochastic process. None the less, to avoid potential bias in the astrometric parameters, direct modelling of the solar-wind effect should be recommended in the timing analysis, unless the anticipated scale of the effect is clearly well below the sensitivity of the data set (see Section 5 for more discussions).

### 3.3 Fitting for $\bar{\omega}$ and $n_{sw}$ in the timing data

Here we conduct additional simulations to investigate the measurability of $\bar{\omega}$ and $n_{sw}$. In these simulations, both parameters are fitted for in the timing analysis, employing the delay model in equations (5) and (6). Since $\Delta_p$ and $\Delta_{sol}$ have a periodic feature with a constant amplitude, the fractional measurement uncertainties of $\bar{\omega}$, $n_{sw}$ is expected to scale with respect to the timing precision ($\sigma_{TOA}$), length ($T$), cadence ($C$) of the data as (Blandford & Teukolsky 1976)

$$ \frac{\sigma_{\bar{\omega}}}{\bar{\omega}} \propto \sigma_{TOA} \times \sqrt{\frac{1}{TC}} \times \frac{1}{\bar{\omega}}, \qquad (11) $$

$$ \frac{\sigma_{n_{sw}}}{n_{sw}} \propto \sigma_{TOA} \times \sqrt{\frac{1}{TC}} \times \frac{1}{n_{sw}}. \qquad (12) $$

Accordingly, for the simulations in this subsection below, we assume 200 ns timing precision, 5 yr of observations and two TOAs per week without loss of generality. The uncertainties of $\bar{\omega}$ and $n_{sw}$ from the simulations can then be scaled with respect to the relations above, to obtain those employing different observational settings. The chosen observing frequency is specified separately in each individual simulation set.

As discussed above, the solar-wind delay has a peculiar deterministic feature in the timing data so is measurable even with single-frequency data. Fig. 9 shows the fractional measurement uncertain-

ties of $\bar{\omega}$ and $n_{sw}$ as a function of ecliptic latitude, each for three choices of typical values. Clearly, the measurement uncertainties of both parameters varies strongly with ecliptic latitude and are lower for small $\beta$ values. Notably, in cases where $n_{sw}$ is significantly measured, i.e. with a fractional error less than 20 per cent, parallax is almost always well measured as well. For $\beta \gtrsim 30$ deg $n_{sw}$ is likely not well constrained, while it is still possible to measure $\bar{\omega}$ with a good precision. This is line with the finding by Tiburzi et al. (2019) that it would be difficult for observations at canonical timing frequencies ($\sim 1$ GHz) to measure solar-wind effect for pulsars at moderate and ecliptic latitudes.

For better demonstration of the fitting result when $\bar{\omega}$ is well measured and $n_{sw}$ is not, Fig. 10 shows the joint posterior distribution of $\bar{\omega}$ and $n_{sw}$ obtained from a Bayesian timing analysis using the TEMPONEST software package. Here, two cases were considered: one $\beta = 3$ deg, $\bar{\omega} = 0.2$ mas, and the other with $\beta = 45$ deg, $\bar{\omega} = 5$ mas, respectively. The first case serves as a reference where both $n_{sw}$ and $\bar{\omega}$ are well measured. In the second case, $n_{sw}$ is not measured significantly and the constraint of its value placed by the Bayesian timing analysis more depends on the chosen prior range. None the less, the global fit including both $n_{sw}$ and $\bar{\omega}$ can still deliver a precise measurement of $\bar{\omega}$, even with a conservative choice of prior range for $n_{sw}$ (see Fig. 10). Using a more informative (and astrophysical) prior slightly improves the measurement precision of $\bar{\omega}$. Although in this scenario applying a fixed $n_{sw}$ in the timing analysis may become feasible, as it may introduce only a minimal bias in the measured parallax as shown in e.g. Fig. 6, conducting the joint fit analysis would still provide useful information for deciding the value of the fixed $n_{sw}$.

The above studies show that for favourable timing precision and ecliptic latitude of the pulsar, the astrometric and solar-wind parameters can be measured jointly even if the data are collected at a single frequency. This means that for data set with poor multifrequency coverage (which could be the case for legacy pulsar timing data set), fitting for a constant solar-wind density could still be attempted so as to avoid the bias in the estimate of astrometric parameters. Of course, collecting data from multifrequency observations (usually the case in modern pulsar timing experiments) could mitigate the correlation between $\bar{\omega}$ and $n_{sw}$ and improve their measurement precision. For a detailed investigation along this line, we simulated data sets each containing two frequency bands. In each iteration, one data set were generated centred at 1.4 GHz and the other at an alternative frequency (which varies between iterations). Fig. 11 shows the measurement uncertainties of $\bar{\omega}$ from a timing fit to the data for various ecliptic latitudes of the pulsar and the frequency ratios of the two bands. The results demonstrate a clear increase in the measurement precision of $\bar{\omega}$ as the frequency separation between the two bands becomes larger, compared to the case where the same amount of observations are all spent at a single frequency. The extent of the improvement depends on the ecliptic latitude of the pulsar. For lower ecliptic latitude, e.g. $\beta < 10$ deg, the precision improves marginally with the extension of frequency coverage, up to a level of approximately 20 per cent. The improvement is more prominent for higher ecliptic latitude cases. For instance, at $\beta$ of 30 deg, the decrease in uncertainty can exceed 40 per cent, while at $\beta$ of 60 deg, it can be as much as 60 per cent when the frequency separation of the two bands is only 20 per cent. This demonstrates that conducting multifrequency observations is more important for pulsars at high ecliptic latitudes, so as to yield a better precision in parallax measurement.

Fig. 12 illustrates the improvement in constraining $\bar{\omega}$ and $n_{sw}$ as the frequency coverage is increased, by showing the joint posterior $1\sigma$ contours of $\bar{\omega}$ and $n_{sw}$ from the Bayesian timing analysis. It can be seen that the contours become tighter and more orthogonal
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Posterior distribution of $\bar{\omega}$ and $n_{sw}$ measured from Bayesian timing analysis of simulated data. The red stars represent the used $n_{sw}$ and $\bar{\omega}$ to simulate the data; these are $\beta = 3$ deg, $\bar{\omega} = 0.2$ mas, $n_{sw} = 8 \text{ cm}^{-3}$ for the upper panel, and $\beta = 45$ deg, $\bar{\omega} = 5$ mas, $n_{sw} = 8 \text{ cm}^{-3}$ for the middle and lower panel. For the analysis in the upper and middle panels, the prior range of $n_{sw}$ was chosen to be [$\bar{n}_{sw} - 20\sigma_0, \bar{n}_{sw} + 20\sigma_0$], where $\bar{n}_{sw}$ is the measurement value from the least-square fit and $\sigma_0$ is its uncertainty. In the case of the middle panel, the prior starts from well below zero to nearly -200 cm$^{-3}$. The prior range used in the analysis of the lower panel is [0, 30] cm$^{-3}$, well including all-known measurements of $n_{sw}$ by far. The ranges of the X and Y axes in the middle and lower panels were set to the same for comparison purpose. The solid, dashed, and dotted lines represent $1\sigma$, $2\sigma$, $3\sigma$ contours, respectively.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 11.** Heat map showing the measurement uncertainties of $\bar{\omega}$ from simulated data consisting of two frequency bands. The simulations were conducted for different ecliptic latitudes and ratios of the frequency of the two bands. The errors were normalized by that at frequency ratio of 1 and with the same $\beta$ angle. Here, we assumed $\bar{\omega} = 1$ mas and $n_{sw} = 8 \text{ cm}^{-3}$. The dotted horizontal line in the lower panel represents frequency ratio of unity and thus normalized uncertainty of 1.</figcaption></figure>

as the frequency ratio of the two bands deviates from unity. This is a result of the decrease in the correlation between these two parameters, as expected from Fig. 4. The reduced correlation allows for more precise and independent measurements of $\bar{\omega}$ and $n_{sw}$, highlighting the benefits of increased frequency separation in multifrequency observations. It is also notable that the improvement on the parameter constraint is more significant for $\beta = 30$ deg than $\beta = 3$ deg. This echoes the conclusion above that multifrequency observations are crucial for precisely measuring $\bar{\omega}$ and $n_{sw}$, if the pulsar is at a high ecliptic latitude.

The above simulations assume equal timing precision at the two frequency bands, while in reality it is likely to be different. Fig. 13 shows the results of the simulation where the alternative band has a different rms residual. Clearly, if the timing precision at the alternative band is higher, there is always improvement in the measurement precision of the parallax, compared to the case of having the same frequency for the two bands. However, if the timing precision of the alternative band is lower, for pulsars at a moderate ecliptic latitude (here $\beta = 45$ deg) or higher, the parallax measurement could still be significantly improved. Take the reference frequency of 1.4 GHz as an example from the figure. For a higher frequency of 2.5 GHz, the uncertainty in $\bar{\omega}$ is still lower even if the timing precision is a factor of a few worse. For frequencies lower than 900 MHz, the precision in $\bar{\omega}$ is higher for rms up to an order of magnitude higher.

## 3.4 A note on observations of pulsars at a very low ecliptic latitude

The analysis above assumes a complete sampling of the solar-wind delay feature around the solar conjunction, where the delay signal is the most prominent. This is a realistic assumption except for pulsars at very low ecliptic latitudes, where practical limitations arise due to the minimum possible solar angle of the source during observations. First, the angular size of the Sun from Earth is 0.524–
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Comparison of joint posterior $1\sigma$ contours of $\bar{\sigma}$ and $n_{sw}$ for data sets with different frequency ratio of the two frequency bands. The upper and lower panels show the results for ecliptic latitude of 3 and 30 deg, respectively. The stars represent the injected value of the two parameters, $\bar{\sigma} = 1$ mas and $n_{sw} = 8$.</figcaption></figure>

0.542 deg. This means that for pulsars at $\beta < 0.3$ deg, they cannot be observed for a short period of time (~2 d) across the solar conjunction because they will be behind the Sun. Additionally, reflecting radio telescopes often have a limit on how close they can point to the Sun to avoid overheating in the receiver focal cabin due to the concentration of optical and infrared radiation.<fnref n="1" /> This constraint means data cannot be collected when the pulsar line of sight is too

<figure><img src="image_18.png" /><figcaption>**Figure 13.** Heat map showing the normalized measurement uncertainties of $\bar{\sigma}$ from simulated data consisting of two frequency bands, with different frequencies and timing precision. The uncertainties of $\bar{\sigma}$ were relative to that at frequency and rms ratio of unity. Here, we assumed $\bar{\sigma} = 1$ mas and $n_{sw} = 8 \text{ cm}^{-3}$, $\beta = 45$ deg. The dotted horizontal and vertical lines represent frequency and rms ratio of unity, respectively.</figcaption></figure>

close to the Sun. Some telescopes that are only able to reflect and concentrate long-wavelength electromagnetic signals in the radio spectrum but not optical and infrared light, would be able to observe closer to the Sun. However, the limiting factor in this case is that the radio signal from the Sun may be included in the side-lobe of the beam, significantly contributing to system temperature and reducing observation sensitivity. For instance, the LOw Frequency ARray (LOFAR) is able to collect useful data for angular separation of the pulsar to the Sun down to approximately 1 deg (Tiburzi et al. 2021).

The lack of data samples around solar conjunction likely impacts the correlation between astrometric and solar-wind parameters and thus affects their measurement precision during practical data analysis. Fig. 14 shows the correlation coefficient of $\bar{\sigma}$ and $n_{sw}$, obtained from data simulated with a given solar-angle limit and at various small ecliptic latitudes. Clearly, in all cases the correlation begins to increase when the solar-angle limit is larger than the ecliptic latitude of the pulsar, resulting in incomplete sampling of the solar-wind delay feature. This leads to less precise measurements of both astrometric and solar-wind parameters compared to the idealistic scenario. It is also worth noticing that for a certain solar-angle limit, $\rho_{\bar{\sigma},n_{sw}}$ is approximately the same for all ecliptic latitudes smaller than this limit. This is likely because the solar-wind delay waveforms in these case are all very similar, once features around the solar conjunction are excluded.

## 4 FITTING FOR SOLAR-WIND DENSITY WITH THE EPTA DR2 DATASET

Following the findings from simulation studies presented above, in this section we reprocessed data from the second data release of the EPTA, to incorporate modelling of $n_{sw}$ in the timing analysis. We used the version ‘DR2full’ data set (see appendix A of EPTA Collaboration 2023a) that contains data from both the legacy and

<footnote n="1">¹For instance, the limit is approximately 3 deg for observations at 1.4 GHz with the Effelsberg Radio Telescope.</footnote>
<figure><img src="image_19.png" /><figcaption>**Figure 14.** Correlation coefficient of $\bar{\omega}$ and $n_{sw}$ as a function of telescope solar angle constraint, for sources at various ecliptic latitudes.</figcaption></figure>

new-generation pulsar backends in the EPTA, mostly at gigahertz frequencies. The analysis focused on a selection of pulsars where $\bar{\omega}$ was measured and the anticipated scale of the solar-wind delay is comparable to the achieved timing precision, allowing for reasonable constraints on $n_{sw}$. The deterministic timing model described in EPTA Collaboration (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 (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

$$S(\ell) = \frac{A^2}{C_0} \left( \frac{\ell}{\ell_t} \right)^{-\gamma}, \qquad (13)$$

where $S(\ell)$, $A$, $\gamma$, $\ell_t$ are the power spectral density as a function of Fourier frequency $\ell$, spectral amplitude, spectral index, and the reference Fourier frequency (set to $1 \text{ 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_0$ equals to $12\pi^2$ for scattering variation and achromatic red noise, and is $k_{\text{DM}}^2$ for DM variation where $k_{\text{DM}} = 2.41 \times 10^{-4} \text{ cm}^{-3} \text{ pc MHz}^2 \text{ 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 $\bar{\omega}$ and $n_{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 $\bar{\omega}$ and $n_{sw}$ presented in Fig. 15.

In the analysis of EPTA Collaboration (2023a), $n_{sw}$ was fixed to $7.9 \text{ cm}^{-3}$ for all pulsars (Madison et al. 2019), except for PSR J1022 + 1001 where a fitted value of $11.1 \pm 0.3 \text{ cm}^{-3}$ was obtained. For PSRs J0030+0451, J0751 + 1807, J1730 - 2304, our timing analysis yields values very close to the previously fixed $n_{sw}$. The measured positions and parallaxes are highly consistent with

those from EPTA Collaboration (2023a), though of slightly lower precision. This outcome is entirely expected due to the inclusion of $n_{sw}$ in the timing fit and the resulting correlation between $\bar{\omega}$, $n_{sw}$.

Meanwhile, for PSRs J1600 - 3053, J1713 + 0747, J1909 - 3744 and J1744 - 1134, the fitted $\bar{\omega}$ and $\beta$ values from our analysis are smaller than those obtained in EPTA Collaboration et al. (2023a), and their fitted $n_{sw}$ values are also correspondingly less than $7.9 \text{ cm}^{-3}$ as anticipated from our simulations. In fact, the inclusion of $n_{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_{sw}$ of $4 \text{ cm}^{-3}$ was used in the timing analysis. For instance, in the case of PSR J1713 + 0747, the fitted value of $n_{sw}$ is only marginally constrained at $4.0 \pm 1.3 \text{ cm}^{-3}$. At the same time, the fitted $\bar{\omega}$ 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 Fig. 6 for a $-3.9 \text{ cm}^{-3}$ offset in $n_{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 (2023a).

The relation between our obtained $n_{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' data set show good consistency with those reported by Susarla et al. (2024) using low-frequency ($\sim 200$ MHz) observations with LOFAR. None the less, 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 (2023b). Notably, the achromatic red-noise measurements remain unchanged after the inclusion of $n_{sw}$ in the timing fit. Therefore, we do not expect significant impact on the gravitational-wave search results presented in EPTA Collaboration (2023c) by fixing the $n_{sw}$ value.

# 5 CONCLUSIONS AND DISCUSSIONS

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_{sw}$. This correlation strongly depends on the ecliptic latitude of the pulsar. We showed that using a fixed $n_{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 data set. 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
<table><caption>**Table 1.** Measurement of timing and noise parameters for a selection of pulsars from our analysis and EPTA Collaboration (2023a) (with subscript 'DR2'), with solar-wind density as a deterministic timing parameter. From the top to bottom, the raws are ecliptic latitude ($\beta$, $\beta_{DR2}$), parallax ($\omega$, $\omega_{DR2}$), solar-wind density ($n_{sw}$), amplitude and spectral index of DM variation ($A_{DM}$, $\gamma_{DM}$), achromatic red noise ($A_{red}$, $\gamma_{red}$), and scattering variation ($A_{Scat}$, $\gamma_{Scat}$). The parallax measurements from EPTA Collaboration (2023a) ($\omega_{DR2}$) were also presented as a comparison. The values in the table represent the median of the posterior distribution from the Bayesian timing analysis with $1\sigma$ credible interval.</caption><thead><tr><th>PSR Jname</th><th>J0030 + 0451</th><th>J0751 + 1807</th><th>J1600 - 3053</th><th>J1713 + 0747</th><th>J1730 - 2304</th><th>J1744 - 1134</th><th>J1909 - 3744</th></tr></thead><tbody><tr><td>β (deg)</td><td>1.4457(14)</td><td>-2.8075798(2)</td><td>-10.07183448(3)</td><td>30.700361492(5)</td><td>0.189(14)</td><td>11.80520318(3)</td><td>-15.15550397(1)</td></tr><tr><td>β<sub>DR2</sub> (deg)</td><td>1.4457(14)</td><td>-2.8075795(2)</td><td>-10.07183439(4)</td><td>30.700361485(5)</td><td>0.189(14)</td><td>11.80520311(3)</td><td>-15.15550393(1)</td></tr><tr><td>ω (mas)</td><td>3.13(7)</td><td>0.68(5)</td><td>0.57(3)</td><td>0.83(2)</td><td>2.09(7)</td><td>2.43(3)</td><td>0.86(2)</td></tr><tr><td>ω<sub>DR2</sub> (mas)</td><td>3.09(6)</td><td>0.85(4)</td><td>0.72(2)</td><td>0.88(1)</td><td>2.08(6)</td><td>2.58(3)</td><td>0.94(2)</td></tr><tr><td>n<sub>sw</sub></td><td>8.4(4)</td><td>5.0(4)</td><td>2.9(4)</td><td>4.0(1.3)</td><td>8.1(7)</td><td>3.2(5)</td><td>3.2(3)</td></tr><tr><td>A<sub>DM</sub></td><td>-</td><td>-11.8<sup>+0.1</sup><sub>-0.1</sub></td><td>-12.2<sup>+0.2</sup><sub>-0.2</sub></td><td>-11.88<sup>+0.04</sup><sub>-0.04</sub></td><td>-11.6<sup>+0.3</sup><sub>-0.3</sub></td><td>-11.84<sup>+0.04</sup><sub>-0.04</sub></td><td>-12.12<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>γ<sub>DM</sub></td><td>-</td><td>3.0<sup>+0.3</sup><sub>-0.3</sub></td><td>4.1<sup>+0.5</sup><sub>-0.5</sub></td><td>1.67<sup>+0.02</sup><sub>-0.02</sub></td><td>2.2<sup>+0.8</sup><sub>-0.8</sub></td><td>1.34<sup>+0.16</sup><sub>-0.16</sub></td><td>1.85<sup>+0.16</sup><sub>-0.15</sub></td></tr><tr><td>A<sub>red</sub></td><td>-14.9<sup>+0.5</sup><sub>-0.6</sub></td><td>-</td><td>-</td><td>-14.2<sup>+0.2</sup><sub>-0.2</sub></td><td>-</td><td>-15.4<sup>+1.2</sup><sub>-1.1</sub></td><td>-15.0<sup>+0.4</sup><sub>-0.5</sub></td></tr><tr><td>γ<sub>red</sub></td><td>5.42<sup>+1.00</sup><sub>-0.87</sub></td><td>-</td><td>-</td><td>3.24<sup>+0.57</sup><sub>-0.52</sub></td><td>-</td><td>5.4<sup>+1.7</sup><sub>-1.8</sub></td><td>5.2<sup>+1.1</sup><sub>-0.9</sub></td></tr><tr><td>A<sub>Scat</sub></td><td>-</td><td>-</td><td>-13.36<sup>+0.04</sup><sub>-0.04</sub></td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>γ<sub>Scat</sub></td><td>-</td><td>-</td><td>1.40<sup>+0.14</sup><sub>-0.14</sub></td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

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 multifrequency 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' data set 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_{sw}$ differed from the fixed value of $7.9 \text{ 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_{sw}$ measurements align well with theirs.

As pointed out by various previous work, the spherical model in equation (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-yr) 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 equation (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. None the less, 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 $\omega$ and $n_{sw}$ strongly depends on several variables, including the pulsar's ecliptic latitude, observing frequency, and timing precision of the data. 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 µ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 $\omega$ is well measured with the data but $n_{sw}$ may not. In this scenario, while fitting for $n_{sw}$ still remains a conservative option, using a given $n_{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 \text{ cm}^{-3}$) lower than the gross average of many PTA pulsars ($7.9 \text{ cm}^{-3}$) obtained in Madison et al. (2019), would likely provide a more realistic representation of the actual $n_{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 modelling 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
<figure><img src="image_20.png" /><figcaption>**Figure 15.** Joint posterior distributions of $\omega$ and $n_{sw}$ measured from EPTA 'DR2full' data set. The stars in the panels correspond to the maximum likelihood value of the parameters.</figcaption></figure>

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 time-scales, 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 canon-

ical timing observations around gigahertz radio frequencies with weekly cadence. We emphasize the significant value of ongoing low-frequency (< 800 MHz) pulsar timing programmes, such as those conducted with the Canadian Hydrogen Intensity Mapping Experiment (CHIME) Telescope (CHIME/Pulsar Collaboration 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.
## ACKNOWLEDGEMENTS

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). AC 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).

## DATA AVAILABILITY

The EPTA 'DR2full' data set are available at: https://zenodo.org/records/8300645.

## REFERENCES

* Agazie G. et al., 2023, ApJ, 951, L9

* Babak S. et al., 2016, MNRAS, 455, 1665

* Blandford R., Teukolsky S. A., 1976, ApJ, 205, 580

* Brisken W. F., Benson J. M., Goss W. M., Thorsett S. E., 2002, ApJ, 571, 906

* CHIME/Pulsar Collaboration, 2021, ApJS, 255, 5

* Cordes J. M., Lazio T. J. W., 2002 preprint (arXiv:astro-ph/0207156)

* Damour T., Deruelle N., 1986, Ann. Inst. H. Poincaré (Physique Théorique), 44, 263

* Deller A. T. et al., 2016, ApJ, 828, 8

* Deller A. T. et al., 2019, ApJ, 875, 100

* Ding H. et al., 2020, MNRAS, 498, 3736

* Ding H. et al., 2023, MNRAS, 519, 4982

* Donner J. Y. et al., 2020, A&A, 644, A153

* EPTA Collaboration, 2023a, A&A, 678, A48

* EPTA Collaboration, 2023b, A&A, 678, A49

* EPTA Collaboration, 2023c, A&A, 678, A50

* Edwards R. T., Hobbs G. B., Manchester R. N., 2006, MNRAS, 372, 1549

* Guo Y. J. et al., 2021, A&A, 654, A16

* Han J. L., Manchester R. N., van Straten W., Demorest P., 2018, ApJS, 234, 11

* Hobbs G., Manchester R., Teoh A., Hobbs M., 2004, in Camilo F., Gaensler B. M., eds, Young Neutron Stars and Their Environments, IAU Symposium 218. Astron. Soc. Pac., San Francisco, p. 139

* Issautier K., Meyer-Vernet N., Moncuquet M., Hoang S., 1998, J. Geophys. Res., 103, 1969

* Iwamoto K. et al., 1998, New Astron., 395, 672

* Jankowski F., van Straten W., Keane E. F., Bailes M., Barr E. D., Johnston S., Kerr M., 2018, MNRAS, 473, 4436

* Keith M. J. et al., 2013, MNRAS, 429, 2161

* Kramer M. et al., 2021, Phys. Rev. X, 11, 041050

* Lam M. T., McLaughlin M. A., Cordes J. M., Chatterjee S., Lazio T. J. W., 2018, ApJ, 861, 12

* Lee K. J., Wex N., Kramer M., Stappers B. W., Bassa C. G., Janssen G. H., Karuppusamy R., Smits R., 2011, MNRAS, 414, 3251

* Lee K. J., Bassa C. G., Janssen G. H., Karuppusamy R., Kramer M., Smits R., Stappers B. W., 2012, MNRAS, 423, 2642

* Lee K. J. et al., 2014, MNRAS, 441, 2831

* Lentati L., Alexander P., Hobson M. P., Feroz F., van Haasteren R., Lee K. J., Shannon R. M., 2014, MNRAS, 437, 3004

* Liu K., Verbiest J. P. W., Kramer M., Stappers B. W., van Straten W., Cordes J. M., 2011, MNRAS, 417, 2916

* Liu K. et al., 2020, MNRAS, 499, 2276

* Madison D. R., Chatterjee S., Cordes J. M., 2013, ApJ, 777, 104

* Madison D. R. et al., 2019, ApJ, 872, 150

* Niu I. C. et al., 2024, MNRAS, 528, 3304

* Parthasarathy A. et al., 2021, MNRAS, 502, 407

* Provornikova E., Opher M., Izmodenov V. V., Richardson J. D., Toth G., 2014, ApJ, 794, 29

* Reardon D. J., Coles W. A., Hobbs G., Ord S., Kerr M., Bailes M., Bhat N. D. R., Venkatraman Krishnan V., 2019, MNRAS, 485, 4389

* Reardon D. J. et al., 2021, MNRAS, 507, 2137

* Shaifullah G., Tiburzi C., Zucca P., 2020, Sol. Phys., 295, 136

* Shklovskii I. S., 1970, Soviet Ast., 13, 562

* Sokól J. M., Bzowski M., Tokumaru M., Fujiki K., McComas D. J., 2013, Sol. Phys., 285, 167

* Splaver E. M., Nice D. J., Stairs I. H., Lommen A. N., Backer D. C., 2005, ApJ, 620, 405

* Susarla S. C. et al., 2024, A&A preprint (arXiv:2409.09838)

* Tiburzi C. et al., 2019, MNRAS, 487, 394

* Tiburzi C. et al., 2021, A&A, 647, A84

* Wang J. B. et al., 2017, MNRAS, 469, 425

* Wang S. Q. et al., 2024, ApJ, 964, 6

* Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29

## APPENDIX: EVALUATION OF THE PARAMETER COVARIANCE

In pulsar timing analysis, the time transformation formula, i.e. the timing model, consists of a sequence of time delay terms that transfers the TOAs of pulsar signals at the observatory to those in the pulsar initial frame, as given by equation (1). The correlation between model parameters can be evaluated by calculating the covariance matrix in the least-square fit by following the method presented in Blandford & Teukolsky (1976) and Damour & Deruelle (1986). Here for simplification, when estimating the correlation between $\omega$ and $n_{sw}$, we only consider the Solar-system Roemer delay term and the two terms introduced by annual parallax and solar wind, respectively. These terms include all timing parameters that are expected to exhibit strong correlation with $\omega$ and $n_{sw}$. Using ecliptic coordinate, namely $\beta$ and $\lambda$, to describe the sky position of the pulsar (and neglecting proper motion), the derivative timing formula can be written as (Blandford & Teukolsky 1976):

$$ \mathcal{R} = -\delta N_0 / \nu + P \cdot \delta\omega + Q \cdot \delta n_{sw} + R \cdot \delta\beta + S \cdot \delta\lambda, \quad (A1) $$

where $P$ and $Q$ are given by equation (8), and $R$ and $S$ are written as

$$ R = \frac{\partial \Delta_{R\odot}}{\partial \beta} + \frac{\partial \Delta_{P}}{\partial \beta} + \frac{\partial \Delta_{\text{sol}}}{\partial \beta}, \quad (A2) $$

$$ S = \frac{\partial \Delta_{R\odot}}{\partial \lambda} + \frac{\partial \Delta_{P}}{\partial \lambda} + \frac{\partial \Delta_{\text{sol}}}{\partial \lambda}. \quad (A3) $$

The derivative terms are expressed by

$$ \frac{\partial \Delta_{R\odot}}{\partial \beta} = \frac{|\vec{r}(E)| \cos (A_T - \lambda) \sin \beta}{c}, \quad (A4) $$

$$ \frac{\partial \Delta_{P}}{\partial \beta} = \frac{\omega |\vec{r}(E)|^2}{2c} \cos^2(A_T - \lambda) \sin 2\beta, \quad (A5) $$
$$ \frac{\partial \Delta_{\text{sol}}}{\partial \beta} = \frac{n_{\text{sw}}}{D_0 |\vec{r}(E)| f^2} \cdot \frac{\sin \beta \cos(A_T - \lambda)(\sin \theta - \theta \cos \theta)}{\sin^3 \theta}, \quad (A6) $$

and

$$ \frac{\partial \Delta_{R\odot}}{\partial \lambda} = \frac{|\vec{r}(E)| \cos \beta \sin(\lambda - A_T)}{c}, \quad (A7) $$

$$ \frac{\partial \Delta_p}{\partial \lambda} = \frac{\omega |\vec{r}(E)|^2}{2c} \cos^2 \beta \sin(2\lambda - 2A_T), \quad (A8) $$

$$ \frac{\partial \Delta_{\text{sol}}}{\partial \lambda} = \frac{n_{\text{sw}}}{D_0 |\vec{r}(E)| f^2} \cdot \frac{\sin(\lambda - A_T) \cos \beta (\sin \theta - \theta \cos \theta)}{\sin^3 \theta}. \quad (A9) $$

The Hessian matrix of the fitted parameters ($\delta N_0$, $\delta\omega$, $\delta n_{sw}$, $\delta\beta$, $\delta\lambda$) is then written as

$$ H = \sum_i \begin{pmatrix} \frac{1}{\nu^2} & -\frac{P_i}{\nu} & -\frac{Q_i}{\nu} & -\frac{R_i}{\nu} & -\frac{S_i}{\nu} \\ -\frac{P_i}{\nu} & P_i^2 & P_i Q_i & P_i R_i & P_i S_i \\ -\frac{Q_i}{\nu} & P_i Q_i & Q_i^2 & Q_i R_i & Q_i S_i \\ -\frac{R_i}{\nu} & R_i P_i & R_i Q_i & R_i^2 & R_i S_i \\ -\frac{S_i}{\nu} & S_i P_i & S_i Q_i & S_i R_i & S_i^2 \end{pmatrix}, \quad (A10) $$

where $\sum_i$ denotes summation over all TOAs. Presume the data have sufficient orbital coverage over the Earth orbit, for a given function $\mathcal{F}$ we can have $\sum_i \mathcal{F}_i \simeq N \cdot \bar{\mathcal{F}}$, where $N$ is the number of TOAs and the overline represents average over the course of an Earth's orbit. The covariance matrix of the fitted parameters is then written as

$$ C = \bar{H}^{-1} = \begin{pmatrix} \frac{1}{\nu^2} & -\frac{\bar{P}}{\nu} & -\frac{\bar{Q}}{\nu} & -\frac{\bar{R}}{\nu} & -\frac{\bar{S}}{\nu} \\ -\frac{\bar{P}}{\nu} & \bar{P}^2 & \bar{P} \bar{Q} & \bar{P} \bar{R} & \bar{P} \bar{S} \\ -\frac{\bar{Q}}{\nu} & \bar{P} \bar{Q} & \bar{Q}^2 & \bar{Q} \bar{R} & \bar{Q} \bar{S} \\ -\frac{\bar{R}}{\nu} & \bar{R} \bar{P} & \bar{R} \bar{Q} & \bar{R}^2 & \bar{R} \bar{S} \\ -\frac{\bar{S}}{\nu} & \bar{S} \bar{P} & \bar{S} \bar{Q} & \bar{S} \bar{R} & \bar{S}^2 \end{pmatrix}^{-1}. \quad (A11) $$

Thus, the correlation coefficients of $\delta n_{sw}$ and $\delta\omega$, $\delta\lambda$, $\delta\beta$, respectively, are given by

$$ \rho_{\delta\omega, \delta n_{\text{sw}}} = \frac{C_{23}}{\sqrt{C_{22} C_{33}}}, \quad (A12) $$

$$ \rho_{\delta\lambda, \delta n_{\text{sw}}} = \frac{C_{34}}{\sqrt{C_{33} C_{44}}}, \quad (A13) $$

$$ \rho_{\delta\beta, \delta n_{\text{sw}}} = \frac{C_{35}}{\sqrt{C_{33} C_{55}}}. \quad (A14) $$

If only considering the two delay terms caused by parallax and solar wind, the derivative timing formula is then written as

$$ R = -\delta N_0 / \nu + P \cdot \delta\omega + Q \cdot \delta n_{\text{sw}}. \quad (A15) $$

Accordingly, the covariance matrix in this case is

$$ C' = \bar{H}^{-1} = \begin{pmatrix} \frac{1}{\nu^2} & -\frac{\bar{P}}{\nu} & -\frac{\bar{Q}}{\nu} \\ -\frac{\bar{P}}{\nu} & \bar{P}^2 & \bar{P} \bar{Q} \\ -\frac{\bar{Q}}{\nu} & \bar{P} \bar{Q} & \bar{Q}^2 \end{pmatrix}^{-1}, \quad (A16) $$

from which gives correlation coefficient of $\delta\omega$ and $\delta n_{sw}$, is then

$$ \begin{aligned} \rho'_{\delta\omega, \delta n_{\text{sw}}} &= \frac{C'_{23}}{\sqrt{C'_{22} C'_{33}}} \\ &= \frac{\overline{H'_{21} H'_{13}} - \overline{H'_{23} H'_{11}}}{\sqrt{(\overline{H'_{33} H'_{11}} - \overline{H'_{31} H'_{13}})(\overline{H'_{22} H'_{11}} - \overline{H'_{21} H'_{12}})}} \\ &= \frac{\overline{P} \cdot \overline{Q} - \overline{P} \cdot \overline{Q}}{\sqrt{\overline{P}^2 - (\overline{P})^2 \sqrt{\overline{Q}^2 - (\overline{Q})^2}}} \end{aligned} \quad (A17) $$

This paper has been typeset from a TeX/LA TeX file prepared by the author.