# Periodic interstellar scintillation variations of PSRs J0613–0200 and J0636+5128 associated with the Local Bubble shell

Yulan Liu, Robert A Main, Joris P.W Verbiest, Ziwei Wu, Krishnakumar M

Ambalappat, Jiguang Lu, David J Champion, Ismaël Cognard, Lucas

Guillemot, Kuo Liu, et al.

► To cite this version:

Yulan Liu, Robert A Main, Joris P.W Verbiest, Ziwei Wu, Krishnakumar M Ambalappat, et al. Periodic interstellar scintillation variations of PSRs J0613–0200 and J0636+5128 associated with the Local Bubble shell. Science China Physics, Mechanics & Astronomy, 2023, 66 (11), pp.119512. 10.1007/s11433-023-2182-6 . hal-04175644

HAL Id: hal-04175644

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

Submitted on 15 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

• Article •

SPECIAL TOPIC:

January 2017 Vol. 60 No. 1: 000000

doi: ****

# Periodic interstellar scintillation variations of PSRs J0613–0200 and J0636+5128 associated with the Local Bubble shell

Yulan Liu<fnref n="1" />,<fnref n="2" />*, Robert A. Main<fnref n="3" />, Joris P. W. Verbiest<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Ziwei Wu<fnref n="1" />,<fnref n="2" />, Krishnakumar M. Ambalappat<fnref n="2" />,<fnref n="3" />, Jiguang Lu<fnref n="1" />, David J. Champion<fnref n="3" />, Ismaël Cognard<fnref n="5" />,<fnref n="6" />, Lucas Guillemot<fnref n="5" />,<fnref n="6" />, Kuo Liu<fnref n="3" />, James W. McKee<fnref n="7" />,<fnref n="8" />, Nataliya Porayko<fnref n="3" />,<fnref n="9" />, Golam. M. Shaifullah<fnref n="9" />,<fnref n="10" />,<fnref n="11" />, and Gilles Theureau<fnref n="5" />,<fnref n="6" />

Received; accepted

Annual variations of interstellar scintillation can be modelled to constrain parameters of the ionized interstellar medium. If a pulsar is in a binary system, then investigating the orbital parameters is possible through analysis of the orbital variation of scintillation. In observations carried out from 2011 January to 2020 August by the European Pulsar Timing Array radio telescopes, PSRs J0613–0200 and J0636+5128 show strong annual variations in their scintillation velocity, while the former additionally exhibits an orbital fluctuation. Bayesian theory and Markov-chain-Monte-Carlo methods are used to interpret these periodic variations. We assume a thin and anisotropic scattering screen model, and discuss the mildly and extremely anisotropic scattering cases. PSR J0613–0200 is best described by mildly anisotropic scattering, while PSR J0636+5128 exhibits extremely anisotropic scattering. We measure the distance, velocity and degree of anisotropy of the scattering screen for our two pulsars, finding that scattering screen distances from Earth for PSRs J0613–0200 and J0636+5128 are 316±²⁸ pc and 262±⁹⁶ pc, respectively. The positions of these scattering screens are coincident with the shell of the Local Bubble towards both pulsars. These associations add to the growing evidence of the Local Bubble shell as a dominant region of scattering along many sightlines.

**Pulsar, Interstellar scintillation, Ionized interstellar medium**

**PACS number(s):** 97.60.Gb, 78.70.Ps, 95.85.Bh

**Citation:** Yulan Liu, Robert A. Main, Joris P. W. Verbiest, et al., *Sci. China-Phys. Mech. Astron.* **60**, 000000 (2017), doi: ****

# 1 Introduction

Pulsars are highly magnetized, rotating neutron stars that emit beams of electromagnetic radiation from their magnetic poles [1] and are embedded in an extremely tenuous inter-

stellar medium (ISM), which contains ordinary matter, relativistic charged particles known as cosmic rays, and magnetic fields [2]. After passing through the warm ionised interstellar medium (IISM), the spatially coherent electromagnetic radiation from the pulsar is distorted, which forms an interfer-

<footnote n="1">¹National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China;</footnote>

<footnote n="2">²Fakultät für Physik, Universität Bielefeld, Bielefeld 33501, Germany;</footnote>

<footnote n="3">³Max-Planck-Institut für Radioastronomie, Bonn 53121, Germany;</footnote>

<footnote n="4">⁴Florida Space Institute, University of Central Florida, Orlando FL 32826, USA;</footnote>

<footnote n="5">⁵Station de radioastronomie de Nançay, Observatoire de Paris, Nançay 18330, France;</footnote>

<footnote n="6">⁶Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans, Orléans 45071, France;</footnote>

<footnote n="7">⁷E.A. Milne Centre for Astrophysics, University of Hull, Hull HU6 7RX, UK;</footnote>

<footnote n="8">⁸Centre of Excellence for Data Science, Artificial Intelligence and Modelling, University of Hull, Hull HU6 7RX, UK;</footnote>

<footnote n="9">⁹Dipartimento di Fisica “G. Occhialini”, Università di Milano-Bicocca, Milano 20126, Italy;</footnote>

<footnote n="10">¹⁰NAF - Osservatorio Astronomico di Cagliari, via della Scienza 5, Selargius (CA) 09047, Italy;</footnote>

<footnote n="11">¹¹INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, Milano 20126, Italy</footnote>
ence pattern at the observer's plane. The interference pattern drifting across the line of sight causes a fluctuation in the source's observed flux density with observing frequency and time, which is defined as interstellar scintillation (ISS) [3]. ISS analysis of pulsars, thus, allows us to probe the smallest-scale distribution and inhomogeneities of the IISM [4].

Scintillation analysis typically relies on measuring the two-dimensional image of the source's observed flux density as a function of observing time and observing frequency, called the dynamic spectrum, where the interference maxima in this dynamic spectrum are called scintles. To quantify the average characteristics of scintles, Cordes & Wolszczan (1986) [5] defined two scintillation parameters: the scintillation bandwidth $v_d$ and scintillation timescale $\tau_d$, where $v_d$ is the half-width at a half maximum along the frequency axis and $\tau_d$ is the half-width at 1/e along the time axis in the two-dimensional auto-correlation function of the dynamic spectrum. At centimetre wavelengths, diffractive interstellar scintillation (DISS) appears as pulsar flux density variations in both time and frequency with characteristic scales ~ minutes and ~ MHz for pulsar timing array (PTA) pulsars, respectively [6, 7]. At metre wavelengths, the scintillation bandwidth is typically as small as ~ kHz for nearby pulsars [8].

The effective velocity, $V_{eff}$, of the line of sight relative to the medium (which is a linear combination of the transverse velocities of Earth, the pulsar, and the IISM) has a strong impact on scintillation parameters, particularly the scintillation timescale $\tau_d$. That is to say, variations in scintillation parameters are the signature of the relative transverse motions of Earth, the pulsar, and the scattering medium. Periodic variations in $\tau_d$ and $v_d$, for instance, are caused by periodically changing transverse velocities of the pulsar and the Earth. Thus, long-term scintillation analyses with periodic variations provide us a possibility to investigate the IISM properties and orbital parameters [9, 11, 40].

Among the prominent features of the local ISM that could be investigated with pulsar scintillation, is the local (hot) bubble. This is a large void in the local Galactic neighbourhood likely created by past supernova explosions [12, 13]. Most likely the local bubble is filled with a hot (~ $10^6$K) ionised gas that is too tenuous to noticeably affect pulsar radiation, although the density of this gas is still a matter of debate [14-16]. It is very well documented, however, that the local bubble contains some interstellar clouds [17, 18] and that it has a very well-defined boundary [19] that can be studied through pulsar scintillation [4, 18].

PSRs J0613-0200 and J0636+5128 both are observed to have strong annual variations in the scintillation timescale, which was reported in Liu et al. (2022) [7]. The main purpose of this work is to make use of these annual variations

to investigate the IISM properties, including the position, the velocity and the anisotropy of the scattering region. Additionally, PSRs J0613-0200 and J0636+5128 are both in a binary system. Another purpose in this work is to determine the orbital inclination angle $i$ and the longitude of ascending node $\Omega_{asc}$ in the pulsar binary system because these two parameters are important for constraining neutron star masses [20] but are difficult to measure through pulsar timing alone. The low-mass companion star and short orbit of PSR J0636+5128 only contribute a small amount to $V_{eff}$, which leads to a negligible orbital period fluctuation in the scintillation parameters [21]. Therefore, we only investigate the orbital parameters for PSR J0613-0200.

In this paper, we present the relevant scintillation background and construct our models assuming an anisotropic thin scattering screen in Section 2. In Section 3, we describe the pulsars and data information, and how we estimate the earth and pulsar velocities. We describe the Bayesian inference and the Markov chain Monte Carlo method that are used to fit the parameters in Section 4. The results and discussion are presented in Section 5. Our conclusions are in Section 6.

## 2 Scintillation from a thin screen

Firstly, we assume that the scintillation is caused by a thin screen at distance $D_s$ from the Earth. For a thin screen model, Cordes & Rickett (1998) [22] presented that the scintillation velocity $V_{ISS}$ can be estimated from the scintillation parameters,

$$V_{ISS} = A_{ISS} \frac{\sqrt{Dv_d}}{f\tau_d}, \qquad (1)$$

where $D$ is the pulsar distance in kpc, $f$ is the observing frequency in GHz, $v_d$ is in units of MHz and $\tau_d$ is in units of seconds.

The factor $A_{ISS}$ depends on the assumed geometry of the scattering medium and the exponent $\alpha$ in the phase structure function on the scattering screen. For a thin screen with Kolmogorov turbulence, $\alpha = 5/3$, and $A_{ISS} \approx 2.78 \times 10^4 \sqrt{2D_s/(D - D_s)}$, where $D_s$ is in units of kpc. However, this formulation is true for isotropic scattering, while scattering is often seen to be quite anisotropic [23]. Stinebring et al. (2022) [24] presented that about 20% of the pulsars in their sample are observed with reverse arclets and a deep valley along the delay axis in the secondary spectrum, indicating substantial anisotropy of scattering. In this work, we consider an anisotropic scattering screen and anisotropy dependent $A_{ISS}$. This anisotropy can be determined by the axial ratio $A_r$ assuming the spatial diffraction pattern as an ellipse, where $A_r = 1$ indicates isotropic scattering. According to the physics of the phase structure function for anisotropic
scattering outlined in Rickett et al. (2014) [9], taking a ge-
ometric mean of the spatial scales on the major and mi-
nor axes of anisotropic scattering, we consider that A<sub>ISS</sub> is
2.78 √(A<sub>r</sub> + 1/A<sub>r</sub>)/2 √2D<sub>s</sub>/(D - D<sub>s</sub>) × 10<sup>4</sup>.

In the thin screen model, we can also predict the scintilla-
tion velocity V<sub>ISS</sub> at the observer as

$$
V_{\text{ISS}} = |V_{\text{eff}}| \frac{D}{D - D_s}. \qquad (2)
$$

Thus, Equation 1 and Equation 2 can be restructured as,

$$
\frac{\sqrt{V_d}}{f\tau_d} \equiv Y = \frac{|V_{\text{eff}}|}{A_{\text{ISS}}} \frac{\sqrt{D}}{D - D_s} \qquad (3)
$$

where we combine all the scintillation observables in the def-
inition of Y.

The effective velocity is a combination of the pulsar’s, the
Earth’s and the IISM’s transverse velocities weighted by the
scattering screen distance as,

$$
V_{\text{eff}} = \frac{D - D_s}{D} V_{\text{E}} + \frac{D_s}{D} (V_{\text{P}} + V_{\mu}) - V_{\text{IISM}}, \quad (4)
$$

where $V_E$ is the Earth velocity, $V_P$ is the binary orbital trans-
verse velocity, $V_\mu$ is the pulsar proper motion transverse ve-
locity and $V_{\text{IISM}}$ is the transverse velocity of the IISM [22].

2.1 Mildly anisotropic scattering

In mildly anisotropic scattering, we take the effective
transverse velocity in equatorial coordinates employing the
quadratic form from Rickett et al. (2014) [9],

$$
\begin{align}
|V_{\text{eff}}| &= \sqrt{aV_{\text{eff},\alpha}^2 + bV_{\text{eff},\delta}^2 + cV_{\text{eff},\alpha}V_{\text{eff},\delta}}, \notag \\
V_{\text{eff},\alpha} &= \frac{D-D_s}{D}V_{\text{E},\alpha} + \frac{D_s}{D}(V_{\text{p},\alpha} + V_{\mu,\alpha}) - V_{\text{IISM},\alpha}, \tag{5} \\
V_{\text{eff},\delta} &= \frac{D-D_s}{D}V_{\text{E},\delta} + \frac{D_s}{D}(V_{\text{p},\delta} + V_{\mu,\delta}) - V_{\text{IISM},\delta}, \notag
\end{align}
$$

where $V_{\text{eff},\alpha}$ and $V_{\text{eff},\delta}$ are the components of the effective ve-
locity in right ascension (α) and in declination (δ), respec-
tively. The coefficients a, b, and c are parameterized by the
axial ratio $A_r$ and the orientation angle of the major axis $\psi$
considering the spatial diffraction pattern as an ellipse,

$$
\begin{align}
a &= \frac{[1 - R \cos(2\psi)]}{\sqrt{1 - R^2}}, \notag \\
b &= \frac{[1 + R \cos(2\psi)]}{\sqrt{1 - R^2}}, \tag{6} \\
c &= \frac{-2R \sin(2\psi)}{\sqrt{1 - R^2}}, \notag
\end{align}
$$

where the orientation angle $\psi$ is clockwise from the right
ascension of the pulsar and $R$ is a scaled parameter for $A_r$:
$R = (A_r^2 - 1)/(A_r^2 + 1)$, which is bound between 0 (indicating
that the inhomogeneities in the plasma density are isotropic)
and 1 (indicating that the inhomogeneities in the plasma den-
sity are extremely anisotropic).

2.2 Extremely anisotropic scattering

Since Brisken et al. (2010) [25] measured an extremely
anisotropic distribution of images for PSR B0834+06 at
327 MHz using very long baseline interferometry, 1D scat-
tering screen models are assumed in increasing numbers of
studies[10, 27]. In the 1D scattering screen model, the effec-
tive transverse velocity is only sensitive to the major axis on
the sky,

$$
\begin{align}
|V_{\text{eff}}| &= \left| \sqrt{aV_{\text{eff},\alpha}^2 + bV_{\text{eff},\delta}^2 + cV_{\text{eff},\alpha}V_{\text{eff},\delta}} - V_{\text{IISM}} \right|, \notag \\
V_{\text{eff},\alpha} &= \frac{D - D_s}{D} V_{\text{E},\alpha} + \frac{D_s}{D} (V_{\text{p},\alpha} + V_{\mu,\alpha}), \tag{7} \\
V_{\text{eff},\delta} &= \frac{D - D_s}{D} V_{\text{E},\delta} + \frac{D_s}{D} (V_{\text{p},\delta} + V_{\mu,\delta}), \notag
\end{align}
$$

where $V_{\text{IISM}}$ is the IISM velocity on the major axis, the coef-
ficients a, b, and c are re-written as,

$$
\begin{align}
a &= \sin^2 \psi, \notag \\
b &= \cos^2 \psi, \tag{8} \\
c &= -2 \sin \psi \cos \psi. \notag
\end{align}
$$

Extremely anisotropic scattering indicates a large A<sub>r</sub>, in
this case, the factor A<sub>ISS</sub> cannot be inferred by the functions
of Rickett et al. (2014) [9]. We therefore give a free scaling
index K (instead of a fixed constant in the isotropic scattering
model) to A<sub>ISS</sub>, thus, A<sub>ISS</sub> = K √2D<sub>s</sub>/(D - D<sub>s</sub>) × 10<sup>4</sup>.

3 Pulsar and data information

Equation 3 is the embodiment of our model, where the left
side of the equation can be seen as the measurement of Y, and
the right side is the prediction of Y. This work is a deeper
analysis of Liu et al. (2022) [7] that presented the annual
variations of scintillation parameters of PSRs J0613–0200
and J0636+5128. The scintillation parameters (v<sub>d</sub> and τ<sub>d</sub>) in
that literature are estimated by exploiting the auto-correlation
function (ACF) of dynamic spectra: v<sub>d</sub> is the half-width at a
half maximum along the frequency axis and τ<sub>d</sub> is the half-
width at 1/e along the time axis in the two-dimensional ACF.
All observations are from three telescopes: the Effelsberg
100-m Radio Telescope (EFF), the Lovell Radio Telescope
at the Jodrell Bank Observatory (JBO) and the Nançay radio
telescope (NRT). The median and the 5/95 percentiles of the
measured scintillation parameters and more observation in-
formation are listed in Table 1. More details about the data
analysis are described Liu et al. (2022) [7].

Due to the limited frequency resolution, most observations
of PSR J0613–0200 from NRT and EFF only give upper lim-
its for v<sub>d</sub>. After several experiments, we noticed that the con-
tribution from the variation of v<sub>d</sub> to the fitted parameters is
<figure><img src="image_3.png" /><figcaption>**Figure 1** Measurements of $v_d$ as a function of time for PSR J0636+5128. These measurements are based on the same data presented by Liu et al. (2022) [7] but have been revised in order to correct for skewness in the dynamic spectrum's ACF, caused by a transverse phase gradient in the scattering screen. Only observations until MJD 57600 are used in the present work since the observing cadence significantly decreased after that (see Liu et al. (2022) [7]).</figcaption></figure>

<table><caption>**Table 1** Measured scintillation parameters and observation information</caption><thead><tr><th>Pulsar</th><th>Telescope</th><th>f<sub>c</sub><br>(MHz)</th><th>MJD range</th><th>N<sub>obs</sub></th><th>Δf<sub>ch</sub><br>(MHz)</th><th>BW<br>(MHz)</th><th>t<sub>sub</sub><br>(sec)</th><th>t̄<sub>obs</sub><br>(min)</th><th>v<sub>d</sub><br>(MHz)</th><th>τ<sub>d</sub><br>(min)</th></tr></thead><tbody><tr><td rowspan="3">J0613-0200</td><td>EFF</td><td>1349</td><td>55661-58019</td><td>55</td><td>1.56</td><td>140.6</td><td>10</td><td>34</td><td>1.7<sup>3.3</sup><sub>1.4</sub></td><td>11<sup>24</sup><sub>7</sub></td></tr><tr><td>JBO</td><td>1556</td><td>56167-58857</td><td>89</td><td>1</td><td>352.0</td><td>65</td><td>35</td><td>2.3<sup>4.2</sup><sub>1.3</sub></td><td>13<sup>27</sup><sub>8</sub></td></tr><tr><td>NRT</td><td>1484</td><td>55817-58852</td><td>194</td><td>4</td><td>512.0</td><td>61</td><td>50</td><td>4.2<sup>9.0</sup><sub>3.0</sub></td><td>16<sup>40</sup><sub>9</sub></td></tr><tr><td rowspan="3">J0636+5128</td><td>EFF</td><td>1347</td><td>56669-59021</td><td>23</td><td>1.56</td><td>200.0</td><td>10</td><td>28</td><td>6.4<sup>37.5</sup><sub>3.4</sub></td><td>8.8<sup>48.8</sup><sub>4.6</sub></td></tr><tr><td>JBO</td><td>1532</td><td>56655-58120</td><td>35</td><td>2</td><td>384.0</td><td>60</td><td>29</td><td>9.7<sup>17.1</sup><sub>4.8</sub></td><td>8.0<sup>32.6</sup><sub>5.1</sub></td></tr><tr><td>NRT</td><td>1544</td><td>56657-58419</td><td>97</td><td>4</td><td>392.0</td><td>61</td><td>60</td><td>12<sup>40</sup><sub>5</sub></td><td>11<sup>56</sup><sub>5</sub></td></tr></tbody></table>

<footnote n="1">1) $f_c$ is the effective center frequency, $N_{obs}$ is the number of observations, $\Delta f_{ch}$ is the channel bandwidth, BW is the effective bandwidth of the observation after the removal of RFI, $t_{sub}$ is the subintegration length, $\bar{t}_{obs}$ is the mean observation length. $v_d$ and $\tau_d$ are the scintillation bandwidth and timescale, respectively. Values in red indicate that measurements are to be considered as upper (for $v_d$) or lower (for $\tau_d$) limits in Liu et al. (2022) [7].</footnote>
<table><caption>**Table 2** Parameters of MSPs J0613–0200 and J0636+5128</caption><thead><tr><th colspan="3">PSR J0613-0200</th></tr><tr><th>Parameter Name</th><th>Parameter value</th><th>Reference paper</th></tr></thead><tbody><tr><td>Right ascension, α (J2000)</td><td>06:13:43.975672</td><td>Desvignes et al. (2016) [28]</td></tr><tr><td>Declination, δ (J2000)</td><td>-02:00:47.22533</td><td>Desvignes et al. (2016) [28]</td></tr><tr><td>Reference epoch for α and δ,</td><td>55000.0</td><td>Desvignes et al. (2016) [28]</td></tr><tr><td>Proper motion in α, μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>1.822(8)</td><td>Desvignes et al. (2016) [28]</td></tr><tr><td>Proper motion in δ, μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>-10.36(2)</td><td>Desvignes et al. (2016) [28]</td></tr><tr><td>Pulsar distance, D (kpc)</td><td>1.1(2)</td><td>Matthews et al. (2016) [29]</td></tr><tr><td>Orbital period, P<sub>b</sub> (days)</td><td>1.198512575184</td><td>Desvignes et al. (2016) [28]</td></tr><tr><td>Epoch of periastron passage, T<sub>0</sub> (MJD)</td><td>53113.953</td><td>Desvignes et al. (2016) [28]</td></tr><tr><td>Projected semi-major axis, x (s)</td><td>1.09144409</td><td>Desvignes et al. (2016) [28]</td></tr><tr><td>Longitude of periastron, ω (°)</td><td>47.2</td><td>Desvignes et al. (2016) [28]</td></tr><tr><td>Orbital eccentricity, e<sub>b</sub> (10<sup>-6</sup>)</td><td>5.4</td><td>Desvignes et al. (2016) [28]</td></tr><tr><td>Orbital inclination, i (°)</td><td>66<sup>+7</sup><sub>-10</sub></td><td>Fonseca et al. (2016) [30]</td></tr><tr><td>Sine of inclination, sin(i)</td><td>0.94(2)</td><td>Perera et al. (2019) [31]</td></tr></tbody></table><table><thead><tr><th colspan="3">PSR J0636+5128</th></tr><tr><th>Parameter Name</th><th>Parameter value</th><th>Reference paper</th></tr></thead><tbody><tr><td>Right ascension, α (J2000)</td><td>06:36:04.84645</td><td>Stovall et al. (2014) [32]</td></tr><tr><td>Declination, δ (J2000)</td><td>51:28:59.9625</td><td>Stovall et al. (2014) [32]</td></tr><tr><td>Reference epoch for α and δ</td><td>56307.0</td><td>Stovall et al. (2014) [32]</td></tr><tr><td>Proper motion in α, μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>4.3(9)</td><td>Stovall et al. (2014) [32]</td></tr><tr><td>Proper motion in δ, μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>2(1)</td><td>Stovall et al. (2014) [32]</td></tr><tr><td>Pulsar distance, D (kpc)</td><td>1.1±0.25</td><td>Kaplan et al. (2018) [21]</td></tr><tr><td>Orbital period, P<sub>b</sub> (days)</td><td>0.066551340763</td><td>Alam et al. (2021) [33]</td></tr><tr><td>Projected semi-major axis, x (s)</td><td>0.00898636</td><td>Alam et al. (2021) [33]</td></tr><tr><td>Longitude of periastron, ω (°)</td><td>5</td><td>Alam et al. (2021) [33]</td></tr><tr><td>Orbital eccentricity, e<sub>b</sub> (10<sup>-5</sup>)</td><td>1.7</td><td>Alam et al. (2021) [33]</td></tr><tr><td>Orbital inclination, i (°)</td><td>24.3±3.5</td><td>Kaplan et al. (2018) [21]</td></tr></tbody></table>
negligible, whereas the contribution from the average value is large. We therefore use the average value of the $v_d$ of JBO observations, which have the highest frequency resolution and have reliable measures of $v_d$, instead of the instantaneous values for PSR J0613-0200. As the ACF of some observations of PSR J0636+5128 shows a skewness resulting from a transverse phase gradient of the scattering screen, the scintillation bandwidths are somewhat biased [34]. We use Equation A6 of Rickett et al. (2014) [9] to correct the skewness for such ACFs, selecting the best fit by eye, and then re-evaluate $v_d$. The revised values of $v_d$ are plotted in Figure 1. Additionally, as the strength of scintillation for PSR J0636+5128 changed significantly and the sampling density drops sharply after MJD 57600, we only use the scintillation measurements before MJD 57600.

We estimate the uncertainty of Y using the fitting procedure and the statistical uncertainties of scintillation parameters. As described in Liu et al. (2022) [7], the uncertainty derived from the ACF is heavily underestimated, we thus include the noise parameters EFAC (an additional scaling factor) and EQUAD (an additional variance added in quadrature) on errors of Y in the following MCMC fitting. In the following context, we use factors F and Q to denote the noise parameters EFAC and EQUAD, respectively.

On the right side of Equation 3, $V_{\text{eff}}$ contains four terms of transverse velocity: the pulsar's proper motion, the pulsar's binary motion, the Earth's orbital motion, and the velocity of the scattering screen. To obtain Earth's velocity, we use a Cartesian coordinate system with the Solar System barycenter as the origin employing the CALCEPH software package [35], then project it onto equatorial coordinates. The pulsar transverse velocity can be estimated using the proper motion $\mu$ that was measured with high precision using timing, $V_\mu = 4.74\mu D$. Assuming a non-relativistic binary system, we need five Keplerian parameters: $P_b$, $T_0$, $x$, $\omega$ and $e$ and two additional parameters: $i$ and $\Omega_{\text{asc}}$ to estimate the pulsar orbital velocity. These five Keplerian parameters are usually well measured with timing. Firstly, the pulsar mean orbital velocity $V_0$ is given by $V_0 = 2\pi xc/(\sin iP_b(1-e^2)^{1/2})$, $x$ is the projected semi-major axis in units of seconds, $c$ is the speed of light in vacuum [36]. The orbital velocity is projected along and perpendicular to the line of nodes in the plane of the sky, marked as $V_{p,\parallel}$ and $V_{p,\perp}$,

$$ 
\begin{aligned}
V_{p,\parallel} &= -V_0 (e \sin \omega + \sin \phi) \\
V_{p,\perp} &= V_0 \cos i (\text{ecos}\omega + \cos \phi)
\end{aligned} 
\qquad (9) 
$$

where $\phi$ is the orbital phase from the line of nodes, $\phi = \theta + \omega$, $\theta$ is the true anomaly and can be calculated by $P_b$, $T_0$ and $e$. Then, the transverse orbital velocity in the directions of right

ascension and declination are

$$ \begin{pmatrix} V_{p,\alpha} \\ V_{p,\delta} \end{pmatrix} = \begin{pmatrix} \cos \Omega_{\text{asc}} & \sin \Omega_{\text{asc}} \\ -\sin \Omega_{\text{asc}} & \cos \Omega_{\text{asc}} \end{pmatrix} \begin{pmatrix} V_{p,\parallel} \\ V_{p,\perp} \end{pmatrix}. \qquad (10) $$

In the next two subsections, we introduce the two pulsars we study in this work, and describe their parameters.

## 3.1 PSR J0613-0200

PSR J0613-0200 is a binary millisecond pulsar (MSP) discovered by Lorimer (1995) [37]. This pulsar is included in all current ongoing PTA experiments [38]. Main et al. (2020) [10] observed scintillation arcs in the 2D power spectrum of scintillation using data from the Large European Array for Pulsars and the Effelsberg telescope, finding the arc curvature varies annually (although they ignored the orbital contribution).

All known parameters used in this work have been listed in Table 2, including the position, the proper motion and five Keplerian parameters [28], and the pulsar distance [29]. Fonseca et al. (2016) [30] and Perera et al. (2019) [31] determined this pulsar's binary orbital inclination angle $i$ using timing. Since only the value of $\sin(i)$ can be determined using timing, we investigate the sense of inclination (i.e., whether $i > 90^\circ$ or $i < 90^\circ$) in this work. In addition, there has been no published measurement of the longitude of ascending node $\Omega_{\text{asc}}$, which we also investigate in this work. The position and the five Keplerian parameters of the binary system for this pulsar have been determined with high precision, and we therefore fix them in the estimation of $V_{\text{eff}}$. In contrast, the uncertainties of the distance and proper motion are included as Gaussian priors.

In conclusion, the scintillation velocity model of PSR J0613-0200 includes nine free parameters in the mildly anisotropic scattering case: $D_s$, $V_{\text{IISM},\alpha}$, $V_{\text{IISM},\delta}$, $A_R$, $\psi$, $i$, $\Omega_{\text{asc}}$, $F$, $Q$; eight free parameters in the extremely anisotropic scattering case: $D_s$, $V_{\text{IISM}}$, $K$, $\psi$, $i$, $\Omega_{\text{asc}}$, $F$, $Q$.

## 3.2 PSR J0636+5128

PSR J0636+5128 is also a binary MSP, discovered in the Green Bank Northern Celestial Cap Pulsar Survey [32]. The position, the proper motion and the distance of this pulsar come from Stovall et al. (2014) [32], with details shown in Table 2. Similar to PSR J0613-0200, the uncertainties of the pulsar distance and proper motion are included in the estimation of $V_{\text{eff}}$. For this pulsar, the orbit period is comparable to the scintillation timescale, and the companion is very light. Therefore the pulsar orbital velocity is low, which means that scintillation parameters do not display measurable
fluctuations from the orbital velocity, and the orbital parame-
ter $T_0$ has not yet been published. Thus, we ignore the pulsar
orbital velocity in the estimation of $V_{\text{eff}}$. Finally, the scin-
tillation velocity model of PSR J0636+5128 includes seven
free parameters in the mildly anisotropic scattering case: $D_s$,
$V_{\text{IISM},\alpha}$, $V_{\text{IISM},\delta}$, $A_r$, $\psi$, $F$, $Q$; six free parameters in the ex-
tremely anisotropic scattering case: $D_s$, $V_{\text{IISM}}$, $K$, $\psi$, $F$, $Q$.

4 Bayesian inference and the Markov chain Monte Carlo method

We use Bayesian theory to derive the parameters of the inves-
tigation list because it is capable of performing a robust mean
for parameter estimation, model selection and visualization
of parameter correlations [39, 41]. From Bayes theorem we
have the joint posterior for a set of model parameters $\theta$,

$$
p(\theta|d) = \frac{p(d|\theta)p(\theta)}{Z} \qquad (11)
$$

where $p(\theta)$ is the prior probability distribution for parame-
ters, $p(d|\theta)$ is the likelihood function of the data $d$ given a
set of parameters of the model and $Z$ is the fully marginal-
ized likelihood function (alternatively the evidence, the prior
predictive and the marginal density of the data). $Z$ is a nor-
malising constant for the posterior and is often ignored in a
single model.

To calculate numerical approximations of multi-
dimensional integrals, the Markov chain Monte Carlo
(MCMC) method was proposed in the last century. The
MCMC method comprises a class of algorithms for sampling
from high-dimensional probability distributions. Markov
Chain is a mathematical system that experiences transitions
from one state to another according to certain probabilistic
rules. Hastings (1970) [42] presented Monte Carlo sampling
methods using Markov Chains and their applications. In this
work, we use a stable and well-tested Python implementation
**emcee** [43] of the affine invariant ensemble sampler [44] for
the MCMC method. Employing Bayesian theory and **emcee**,
we can produce and draw posterior probability distributions
for all studied parameters.

We give a uniform prior within a physically-motivated
bound region for each parameter: $D_s \in (0, D)$, $V_{\text{IISM},\alpha}$ and
$V_{\text{IISM},\delta} \in [-100, 100]$ km s⁻¹, $A_r \in [1, 100]$, $\psi \in [0°, 180°)$,
$i \in [0°, 180°)$, $\Omega_{\text{asc}} \in [-180°, 180°)$, and run **emcee** with $10^5$
steps. The first 25% of the steps in the sampling process are
considered to be 'burn-in' stage and are therefore discarded.
In that stage, the chain is not effectively exploring the param-
eter space as the global maximum is not yet found.

5 Results and discussion

We show the best-fit results for PSRs J0613–0200 and
J0636+5128 in Figures 2 and 3, respectively, including the
most likely values with uncertainties (the 1-D and 2-D poste-
rior probability distributions) for all fitted parameters.

We use Equation 3 in the MCMC fitting to separate
the scintillation observables from our scintillation velocity
model. In order to show the results from the MCMC fitting
in the physical sense, we plot the scintillation velocities from
Equations 1 and 2 for PSRs J0613–0200 and J0636+5128 in
Figures 4 and 5, respectively. The $V_{\text{ISS,M}}$ is the scintillation
velocity from Equation 1, $V_{\text{ISS,P}}$ is that from Equation 2. We
plot the differences between the two scintillation velocities
and present the reduced $\chi^2$ in the second panel. In the third
and fourth panels, all scintillation velocities are plotted as a
function of the day of year and orbital phase, respectively.

<figure><img src="image_4.png" /><figcaption>**Figure 4** The scintillation velocities of PSR J0613–0200. $V_{\text{ISS,M}}$ is the scintillation velocity from Equation 1, $V_{\text{ISS,P}}$ is that from Equation 2. The $V_{\text{ISS,P}}$ values in the first panel are determined using a cadence of 0.1 days. The second panel shows the differences between $V_{\text{ISS,M}}$ and $V_{\text{ISS,P}}$. In the third and fourth panels, we plot all points as a function of day of year and orbital phase, respectively.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 2** The posterior probability distributions of all fitted parameters for PSR J0613-0200. In the 1-D histograms, the red lines are the kernel-density estimate smoothed versions of the distributions, the vertical dashed lines indicate the 10% fractional percentiles, the most likely values and the 90% fractional percentiles. The most likely values and the upper/lower errors are presented at the top of the 1-D histograms. Factors *F* and *Q* are the noise parameters EFAC and EQUAD, respectively. The factor *Q* is scaled up by a factor of one thousand, and *D*<sub>s</sub> is in units of pc.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3** The posterior probability distributions of all fitted parameters for PSR J0636+5128. In the 1-D histograms, the red lines are the kernel-density estimate smoothed versions of the distributions, the vertical dashed lines indicate the 10% fractional percentiles, the most likely values and the 90% fractional percentiles. The most likely values and the upper/lower errors are presented at the top of the 1-D histograms. Factors *F* and *Q* are the noise parameters EFAC and EQUAD, respectively. The factor *Q* is scaled up by a factor of one thousand, and *D*<sub>s</sub> is in units of pc.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5** The scintillation velocities of PSR J0636+5128. $V_{\text{ISS,M}}$ is the scintillation velocity from Equation 1, $V_{\text{ISS,P}}$ is that from Equation 2. The $V_{\text{ISS,P}}$ values in the first panel are determined using a cadence of 0.1 days. The second panel shows the differences between $V_{\text{ISS,M}}$ and $V_{\text{ISS,P}}$. In the third and fourth panels, we plot all points as a function of day of year and orbital phase, respectively. Since the absence of the epoch of periastron passage $T_0$, $V_{\text{ISS,P}}$ is omitted from the fourth panel.</figcaption></figure>

For PSR J0613-0200, the best fit is from the mildly anisotropic scattering model. In Figure 2, all fitted parameters converge rapidly and most parameters show a well-behaved Gaussian distribution in the posterior probability distribution. The reduced $\chi^2$ in the mildly anisotropic scattering model is 3.42 which is better than that of the extremely anisotropic scattering model (6.52). The best fit for PSR J0636+5128 is from the extremely anisotropic scattering model and is shown in Figure 3, where the reduced $\chi^2$ of the scintillation velocities is 0.91. Also, all fitted parameters for this pulsar converge rapidly and show a well-behaved Gaussian distribution.

## 5.1 Positions and velocities of the scintillation screens

Based on the posterior probability distributions discussed above, the most likely scintillation screen distances for PSRs J0613-0200 and J0636+5128 are, respectively, $316_{-20}^{+28}$ pc$^1$) and $262_{-38}^{+96}$ pc. Both of these values are much closer to the Earth than to the pulsar; and they are in the vicinity of the Local Bubble boundary. In this section we will investigate this possible association.

In Figure 6, we plot three crosscuts of the solar neighborhood based on the 3D dust extinction map of Lallement et al. (2019) [46], which can be assumed a proxy for the gas density map. The three crosscuts clearly show the Local Bubble cavity. We overlay the projections of the scintillation screen positions (with error bars) for our two pulsars on the three crosscuts, which show that the scintillation screens of PSRs J0613-0200 and PSR J0636+5128 align well with the Local Bubble boundary.

Next, we calculate the distances of the inner local bubble surface on the lines of sight to our two pulsars using the Local Bubble inner surface model from Pelgrims et al. (2020) [47], finding a surface distance of 252 pc towards PSR J0613-0200 and 202 pc towards PSR J0636+5128. Both measured screen distances are slightly beyond this boundary, indicating the scattering screens of our pulsars are likely associated with the Local Bubble shell.

In Figure 7, we plot the gas density distribution along the line of sight of two pulsars employing the 3D dust extinction map of Lallement et al. 2019 [46], and mark the positions of the scattering screens and the inner local bubble surfaces. The positions of the scattering screens coincide with the first peaks of the gas density distributions along the lines of sight to PSRs J0613-0200 and J0636+5128, respectively, which supports the hypothesis that the nearest dense IISM structure is the primary cause of the observed interstellar scintillation, assuming a positive correlation between the neutral and ionized parts of the medium. Using this hypothesis, we estimate the lower limit of the shell thickness for the local bubble to be twice the distance between the scattering screen and the inner local bubble surface. In the directions of PSRs J0613-0200 and J0636+5128, the lower limits on the shell thickness are $128_{-20}^{+28}$ pc and $120_{-38}^{+96}$ pc, respectively.

The transverse IISM velocities for PSRs J0613-0200 and J0636+5128 are around 18 km/s and -2 km/s, respectively. These velocities are well within the range of expectations for the ISM residing on the Local Bubble shell. For PSR J0613-0200, the projection of the transverse IISM velocities onto the major axis of the IISM is roughly 12 km/s, which is in agreement with that of the 2013 event reported in

<footnote n="1">1) We note our scintillation screen distance to PSR J0613-0200 is consistent with the one derived by Main et al. (2020) [10] and Main et al. (2023) [45]</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6** Scintillation screen position on three crosscuts (XY, XZ and YZ planes) of the 3D dust extinction map from Lallement et al. 2019 [46]. The X axis points from the Sun (the circle dot) to the Galactic center, the Y axis points towards *l* = 90° and the Z axis points to the North Galactic pole at Galactic latitude *b* = 90°. The colour scale shows log(*A*′<sub>*v*</sub>) and indicates the gas density, where *A*′<sub>*v*</sub> is the differential extinction in units of magnitudes per parsec.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7** The gas density distribution along the line of sight to two pulsars, from the 3D dust extinction map of Lallement et al. 2019 [46]. The coloured vertical lines denote the estimated position of the scattering screen for each pulsar, while the lighter shades of the vertical lines indicate the uncertainty range of the scintillation screen positions. The dashed black vertical lines indicate positions of the inner local bubble surface on the line of sight to each pulsar.</figcaption></figure>
Main et al. (2020) [10]. Main et al. (2023) [45] presented the IISM velocities towards PSR J0613−0200 in both directions of right ascension and declination. In our estimates, the IISM velocity towards PSR J0613−0200 in right ascension is $17.76_{-0.05}^{+2.19}$ km/s, which is consistent with their result. However, in the direction of declination, the IISM velocity in our estimates is $1.22_{-0.90}^{+1.47}$ km/s, which slightly deviates from their result.

## 5.2 Anisotropic scattering

From our modelling, the scattering screens are both evidently anisotropic. The axial ratio $A_r$ of PSR J0613−0200 is $1.6_{-0.2}^{+0.3}$, which is in agreement with the mildly anisotropic scattering model. The orientation of the major axis $\psi$ is $40_{-6}^{+7}$ degrees, which is consistent with the orientation reported for the 2013 event in Main et al. (2020) [10], but differs from the orientations reported in Main et al. (2023)$^{2}$[45]. For PSR J0636+5128, the periodic variation of scintillation parameters is most accurately explained by the extremely anisotropic scattering model and $\psi$ is $83_{-3}^{+4}$ degrees, indicating that the diffraction pattern is a filament oriented parallel to the direction of declination of this pulsar.

Although the posterior distributions of $V_{\text{IISM},\delta}$ and $A_r$ for PSR J0636+5128 do not converge well in the mildly anisotropic scattering model, the scintillation velocities $V_{\text{ISS,M}}$ and $V_{\text{ISS,P}}$ are well matched with a reduced $\chi^2$ of 1.91 for the differences. Due to the negligible fluctuation of scintillation parameters caused by the pulsar orbital velocity and the absence of the epoch of periastron passage $T_0$, we neglected the pulsar orbital velocity when calculating $V_{\text{eff}}$. As a result, we were unable to unambiguously constrain the scattering screen parameters and hence cannot entirely exclude the mildly anisotropic scattering model. This can be examined with well-resolved scintillation arcs, requiring longer observations and finer frequency resolution than we had in this work.

## 5.3 Orbital angles $i$ and $\Omega_{\text{asc}}$

Fonseca et al. (2016) [30] and Perera et al. (2019) [31] published the sine value of the inclination angle of PSR J0613−0200 as $0.91_{-0.08}^{+0.05}$ and $0.94_{-0.02}^{+0.02}$, respectively. From our scintillation analysis, we can break that ambiguity by fitting both the sine and cosine of $i$. In our modelling, the inclination angle $i$ of PSR J0613−0200 is $101_{-18}^{+8}$ degrees, the $\sin(i)$ is $0.98_{-0.03}^{+0.02}$ which is consistent with the timing result, but it points to $i>90^\circ$. Additionally, we present the $\Omega_{\text{asc}}$ of PSR J0613−0200 as $-48_{-9}^{+36}$ degrees. Unfortunately, the values of $i$ and $\Omega_{\text{asc}}$ in our modelling are different with Main

et al. (2023) [45].

# 6 Conclusions

Assuming an anisotropic thin scattering screen model, we have modelled the annual variations of scintillation parameters using Equation 3, allowing us to constrain the position, velocity and anisotropy of scattering screens of PSRs J0613−0200 and J0636+5128. Comparing our results for PSR J0613−0200 with Main et al. (2020) [10] and Main et al. (2023) [45], our estimates for the distance of the scattering screen are consistent with the results in literature. Both scattering screens in this work are consistent with being associated with the Local Bubble shell. Since the Local Bubble shell appears as a closed surface [48], the electron density fluctuations sharply increase at the shell along the propagation path of pulsar signals. Thus, for most nearby pulsars, the Local Bubble shell plays a substantial role in ISS [4, 49, 50].

Early scintillation studies typically assumed isotropic scattering [51]. Over the past two decades, there has been increasing evidence towards anisotropic and inhomogeneous [52] scattering in many cases. Since the discovery of an extremely anisotropic distribution of images by Brisken et al. (2010) [25], the hypothesis of extremely anisotropic scattering has been frequently used in recent years. This hypothesis only requires velocities along the major axis of anisotropy, which results in fewer input parameters being necessary [26]. Extremely anisotropic scintillation is supported in a few cases, e.g., sources with sharp inverted arclets (PSRs B0834+06, B1508+55 and B0450−18 [53-55]). From our modelling, the scattering screens are both evidently anisotropic, one is even extremely anisotropic.

Additionally, we have presented $\Omega_{\text{asc}}$ and $i$ for PSR J0613−0200. The longitude of ascending node $\Omega_{\text{asc}}$ is presented for the first time. The calculation of the pulsar orbital velocity requires the introduction of both sine and cosine functions of $i$, which enable the determination of $i = 101_{-18}^{+8}$ ($> 90^\circ$), thus resolving the ambiguity which timing was unable to resolve.

We thank Bill Coles for many useful discussions and suggestions. This work is supported by the National Natural Science Foundation of China (Grant No. 12003047), the Major Science and Technology Program of Xinjiang Uygur Autonomous Region (No. 2022A03013-2) and Natural Science Foundation of Xinjiang Uygur Autonomous Region (No. 2022D01D85). JPWV acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the Heisenberg programme (Project No. 433075039). Part of this work is based on observations with the 100-m radio telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg in Germany. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a consolidated grant from the STFC in the UK. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche

<footnote n="2">2) Main et al. (2023) [45] defined that $\psi$ is clockwise from the declination of the pulsar.</footnote>
Scientifique (CNRS). We acknowledge financial support from "Programme National de Cosmologie et Galaxies" (PNCG) of CNRS/INSU, France.

1 Gold, T. 1968, Nature, 218, 731

2 Ferrière, K. M. 2001, Reviews of Modern Physics, 73, 1031

3 Scheuer, P. A. G. 1968, Nature, 218, 920

4 Bhat, N. D. R., Gupta, Y., & Rao, A. P. 1998, ApJ, 500, 262

5 Cordes, J. M. & Wolszczan, A. 1986, ApJ, 307, L27

6 Levin, L., McLaughlin, M. A., Jones, G., et al. 2016, ApJ, 818, 166

7 Liu, Y., Verbiest, J. P. W., Main, R. A., et al. 2022, A&A, 664, A116

8 Wu, Z., Verbiest, J. P. W., Main, R. A., et al. 2022, A&A, 663, A116

9 Rickett, B. J., Coles, W. A., Nava, C. F., et al. 2014, ApJ, 787, 161

10 Main, R. A., Sanidas, S. A., Antoniadis, J., et al. 2020, MNRAS, 499, 1468

11 Mall, G., Main, R. A., Antoniadis, J., et al. 2022, MNRAS, 511, 1104

12 Smith, R. K. & Cox, D. P. 2001, ApJS, 134, 283.

13 Heiles, C. 1998, ApJ, 498, 689. doi:10.1086/305574

14 Welsh, B. Y. & Shelton, R. L. 2009, Ap&SS, 323, 1.

15 Galeazzi, M., Chiao, M., Collier, M. R., et al. 2014, Nature, 512, 171.

16 Snowden, S. L., Chiao, M., Collier, M. R., et al. 2014, ApJ, 791, L14.

17 Linsky, J. L. & Redfield, S. 2014, Ap&SS, 354, 29.

18 Linsky, J. L., Rickett, B. J., & Redfield, S. 2008, ApJ, 675, 413.

19 Lallement, R., Vergely, J.-L., Valette, B., et al. 2014, A&A, 561, A91.

20 Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nature, 467, 1081

21 Kaplan, D. L., Stovall, K., van Kerkwijk, M. H., Fremling, C., & Is-
trate, A. G. 2018, ApJ, 864, 15

22 Cordes, J. M. & Rickett, B. J. 1998, ApJ, 507, 846

23 Coles, W. A., McLaughlin, M. A., Rickett, B. J., Lyne, A. G., & Bhat, N. D. R. 2005, ApJ, 623, 392

24 Stinebring, D. R., Rickett, B. J., Minter, A. H., et al. 2022, ApJ, 941, 34

25 Brisken, W. F., Macquart, J. P., Gao, J. J., et al. 2010, ApJ, 708, 232

26 Walker, M. A., de Bruyn, A. G., & Bignall, H. E. 2009, MNRAS, 397, 447

27 Main, R. A., Bethapudi, S., Marthi, V. R., et al. 2022, arXiv e-prints, arXiv:2212.04839

28 Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341

29 Matthews, A. M., Nice, D. J., Fonseca, E., et al. 2016, ApJ, 818, 92

30 Fonseca, E., Pennucci, T. T., Ellis, J. A., et al. 2016, ApJ, 832, 167

31 Perera, B. B. P., DeCesar, M. E., Demorest, P. B., et al. 2019, MNRAS,

490, 4666

32 Stovall, K., Lynch, R. S., Ransom, S. M., et al. 2014, ApJ, 791, 67

33 Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2021, ApJS, 252, 4

34 Yao, J., Zhu, W., Manchester, R. N., et al. 2021, Nature Astronomy, 5, 788

35 Gastineau, M., Laskar, J., Manche, H., & Fienga, A. 2015, CALCEPH: Planetary ephemeris files access code, Astrophysics Source Code Library, record ascl:1505.001

36 Lorimer, D. R. & Kramer, M. 2004, Handbook of Pulsar Astronomy, Vol. 4

37 Lorimer, D. R. 1995, MNRAS, 274, 300

38 Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267

39 Reardon, D. J., Coles, W. A., Hobbs, G., et al. 2019, MNRAS, 485, 4389

40 Reardon, D. J., Coles, W. A., Bailes, M., et al. 2020, ApJ, 904, 2

41 Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, ApJS, 241, 27

42 Hastings, W. K. 1970, Biometrika, 57, 97

43 Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

44 Goodman, J. & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65

45 Main, R. A., Antoniadis, J., Chen, S., et al. 2023, arXiv:2306.13462. doi:10.48550/arXiv:2306.13462

46 Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135

47 Pelgrims, V., Ferrière, K., Boulanger, F., Lallement, R., Montier, L. 2020, A&A, 636, A17

48 Zucker, C., Goodman, A. A., Alves, J., et al. 2022, Nature, 601, 334

49 Bhat, N. D. R., Ord, S. M., Tremblay, S. E., McSweeney, S. J., & Tingay, S. J. 2016, ApJ, 818, 86

50 Bhat, N. D. R. & Gupta, Y. 2002, ApJ, 567, 342

51 Rickett, B. J. 1977, ARA&A, 15, 479

52 Wu, Z., Coles, W. A., Verbiest, J. P. W., et al. 2023, MNRAS, 520, 5536

53 Hill, A. S., Stinebring, D. R., Asplund, C. T., et al. 2005, ApJ, 619, L171

54 Rickett, B. J., Stinebring, D. R., Zhu, H., & Minter, A. H. 2021, ApJ, 907, 49

55 Sprenger, T., Main, R., Wucknitz, O., Mall, G., & Wu, J. 2022, MNRAS, 515, 6198

56 Bialy, S., Zucker, C., Goodman, A., et al. 2021, ApJ, 919, L5