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

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 interstellar 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 interference 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 $\nu_{\rm d}$ and scintillation timescale $\tau_{\rm d}$, where $\nu_{\rm d}$ is the half-width at a half maximum along the frequency axis and $\tau_{\rm 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 $\sim$ minutes and $\sim$ MHz for pulsar timing array (PTA) pulsars, respectively [6, 7]. At metre wavelengths, the scintillation bandwidth is typically as small as $\sim$ kHz for nearby pulsars [8].

The effective velocity, $V_{\rm{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_{\rm 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_{\rm d}$ and $\nu_{\rm 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, 40, 11].

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 ($\sim 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, 15, 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 [18, 4].

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_{\rm 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_{\rm{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.

The factor $A_{\rm 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_{\rm 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_{\rm ISS}$. This anisotropy can be determined by the axial ratio $A_{\rm r}$ assuming the spatial diffraction pattern as an ellipse, where $A_{\rm 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 geometric mean of the spatial scales on the major and minor axes of anisotropic scattering, we consider that $A_{\rm ISS}$ is 2.78$\sqrt{(A_{r}+1/A_{r})/2}\sqrt{2D_{s}/(D-D_{s}})\times 10^{4}$.

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,

where $V_{\rm E}$ is the Earth velocity, $V_{\rm p}$ is the binary orbital transverse velocity, $V_{\rm\mu}$ is the pulsar proper motion transverse velocity and $V_{\rm IISM}$ is the transverse velocity of the IISM [22].

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 ($\nu_{\rm d}$ and $\tau_{\rm d}$) in that literature are estimated by exploiting the auto-correlation function (ACF) of dynamic spectra: $\nu_{\rm d}$ is the half-width at a half maximum along the frequency axis and $\tau_{\rm d}$ 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 information 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 limits for $\nu_{\rm d}$. After several experiments, we noticed that the contribution from the variation of $\nu_{\rm d}$ to the fitted parameters is negligible, whereas the contribution from the average value is large. We therefore use the average value of the $\nu_{\rm d}$ of JBO observations, which have the highest frequency resolution and have reliable measures of $\nu_{\rm 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 $\nu_{\rm d}$. The revised values of $\nu_{\rm 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_{\rm 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_{\rm\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_{\rm 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{x}c/(\sin\,iP_{\rm 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_{\rm p,\parallel}$ and $V_{\rm p,\perp}$,

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_{\rm b}$, $T_{0}$ and $e$. Then, the transverse orbital velocity in the directions of right ascension and declination are

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

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

where $p(\theta)$ is the prior probability distribution for parameters, $p(d|\theta)$ is the likelihood function of the data $d$ given a set of parameters of the model and $Z$ is the fully marginalized likelihood function (alternatively the evidence, the prior predictive and the marginal density of the data). $Z$ is a normalising 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_{\rm IISM,\alpha}$ and $V_{\rm IISM,\delta}\in[-100,100]$ km s^-1, $A_{r}\in[1,100]$, $\psi\in[0^{\circ},180^{\circ})$, $i\in[0^{\circ},180^{\circ})$, $\Omega_{\rm asc}\in[-180^{\circ},180^{\circ})$, and run emcee with 10⁵ 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 parameter space as the global maximum is not yet found.

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 posterior 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_{\rm{ISS,M}}$ is the scintillation velocity from Equation 1, $V_{\rm{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.

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.

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, 50, 49].

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, 54, 55]). From our modelling, the scattering screens are both evidently anisotropic, one is even extremely anisotropic.

Additionally, we have presented $\Omega_{\rm asc}$ and $i$ for PSR J0613$-$0200. The longitude of ascending node $\Omega_{\rm 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^{+8}_{-18}~{}(\textgreater\,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 Xinjiiang 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 Scientifique (CNRS). We acknowledge financial support from “Programme National de Cosmologie et Galaxies” (PNCG) of CNRS/INSU, France.