Observations of Low-frequency Radio Emission from Millisecond Pulsars and Multipath Propagation in the Interstellar Medium

Observational data presented in this paper were obtained during 2016–2017, when the MWA underwent a major upgrade, which involved the deployment and commissioning of 128 new tiles, effectively expanding on the original array, i.e., the Phase I MWA. Of the newly added tiles, 72 are configured into 2 × 36 hexagonal layouts (each termed a "Hex") to provide a large number of redundant baselines, while the remaining 56 are placed farther out from the original array to provide longer baselines up to ∼6 km (i.e., doubling of the imaging resolution and reducing the confusion limit). The basic system architecture, in terms of ability to correlate or voltage-record, is, however, limited to a maximum of 128 tiles at a given time. An important feature is that the array can be periodically reconfigured either as a compact array, in which the majority of the 128 tiles are located within a central region 300 m across (comprising the two Hex's and the core), or as an extended array that includes 56 outer tiles and some fraction of the core. These new capabilities are formally referred to as the Phase II MWA.

Even though the equivalent collecting area, or the achievable sensitivity, remains essentially the same as far as pulsar observations are concerned, calibration and beam-forming with the compact configuration will naturally be far less susceptible to potential residual errors from ionospheric calibration.¹² Such a compact configuration is also extremely appealing for undertaking large-area pulsar surveys given the significant reduction in the computational cost for realizing beam-forming across (or pixelizing) the full field of view (FOV).

VCS-recorded data can be coherently summed to generate a tied-array (phased-array) beam on the sky by performing a coherent addition of voltage signals from individual tiles (S. M. Ord et al. 2018, in preparation). This maximizes the achievable sensitivity and is thus the preferred mode for pulsar observations with the MWA. However, this procedure involves incorporating polarimetric beam responses of individual tiles (using analytic models for the beams), as well as the cable and geometric delay models and complex gain information (amplitude and phase) for each tile. The calibration can make use of an in-field bright source, or multiple moderately bright sources in the field, or sometimes pointed observations of one of the bright calibrators (e.g., Hydra A, 3C 444, Hercules A) taken prior to pulsar observations. In all these cases, first the visibilities are created using the MWA correlator (Ord et al. 2015), and then for each sub-band a calibration solution (for both amplitude and phase) is iteratively generated for each tile using the Real Time System (RTS; Mitchell et al. 2008). The output from the RTS is a calibration solution for each coarse channel, i.e., 24 solutions per observation for each tile. Between iterations, tiles with poor-quality calibration solutions are flagged out (the number tends to vary depending on the observation). Once solutions have converged for each sub-band, a coherent tied-array beam is produced using only those tiles with good solutions. The full processing chain runs on the Galaxy cluster of the Pawsey supercomputing facility and is routinely employed for sensitive pulsar observations (e.g., McSweeney et al. 2017; Meyers et al. 2017).

There is, however, an important caveat associated with the choice of the in-field calibration strategy. Since the tile (primary) beam pattern is a strong function of the observing frequency, if the calibrator is not in close proximity to the pulsar, it may lead to difficulties calibrating higher-frequency bands as the calibrator becomes near (or outside) the edge of the beam (and hence not ideal for reliable calibration). For example, Pictor A, which is located at ∼10° from the position of PSR J0437−4715, is no longer a suitable in-field calibrator for this pulsar at frequencies ≳220 MHz. Observations across the full MWA band may therefore warrant considering additional (or alternate) strategies to ensure successful calibration across the full frequency range.

The output of the second-stage beam-former in the EDA signal path is low-pass filtered and sampled at a rate of 655.36 MHz. For pulsar observations, the resultant data streams are rescaled to 4-bit samples before writing out to disks as a dual-polarization voltage time series. These data can therefore be processed using the phase-coherent dispersion removal technique over the majority of the 327.68 MHz band, enabling pulsar observations over a large instantaneous bandwidth. In practice, the available bandwidth is limited by the response range of MWA dipoles (i.e., ∼50–300 MHz) and the rapid growth of dispersive delay at low radio frequencies. The DSPSR software package,¹⁵ which has been the backbone of multiple generations of pulsar back ends at the Parkes 64 m telescope, has been integrated with the output of the EDA data stream. The baseband data stream can be coherently de-dispersed while dividing the full EDA band into sub-bands (typically 32k channels) using the convolving filter bank (van Straten & Bailes 2011) to generate virtually artifact-free pulse profiles. For a given DM, the lowest radio frequency amenable to coherent de-dispersion depends on the computational memory available for processing.

Flux calibration for MWA pulsar detections made use of the technique developed by Meyers et al. (2017). In summary, we simulated the tied-array (phased-array) beam pattern by modeling it as the product of the tile beam pattern and the array factor (see their Equations (11) and (12)). The tile beam pattern was simulated using the formalism as per Sutinjo et al. (2015), and the array factor depends on the configuration used for the observation (i.e., whether Phase I or Phase II, and the number of tiles used to form a tied-array beam on the pulsar). The simulated tied-array beam is then convolved with the global sky model (de Oliveira-Costa et al. 2008) and integrated over the sky (see Sokolowski et al. 2015) to produce estimates of the antenna temperature ${T}_{\mathrm{ant}}$. These are then used to calculate the system temperature ${T}_{\mathrm{sys}}$, under the assumption that the radiation efficiency η_rad ≃ 1 and using the receiver temperature measurements (T_rcv ≈ 40 K) for the MWA. The integral of this array factor power pattern is then used to determine the beam solid angle ${{\rm{\Omega }}}_{{\rm{A}}}$, from which the tied-array gain $G={A}_{{\rm{e}}}/2{k}_{{\rm{B}}}$ (where A_e is the effective collecting area and ${k}_{{\rm{B}}}$ is the Boltzmann constant) is estimated and translated into the system equivalent flux density SEFD = ${T}_{\mathrm{sys}}$/G. Using this approach, we estimated G and ${T}_{\mathrm{sys}}$ for each sub-band of observations across the MWA band, by factoring in also the directional and frequency dependencies of the MWA tile pattern. We then use the standard radiometer equation for pulsar observations (i.e., accounting for the pulsar duty cycle W/P, where W is the equivalent pulse width and P is the pulsar period) to calculate the corresponding flux density scales for our observations (Lorimer & Kramer 2004). The equivalent width W is estimated as ${A}_{\mathrm{pulse}}$/${S}_{\mathrm{peak}}$, where ${A}_{\mathrm{pulse}}$ is the integrated flux within the "on" pulse region and ${S}_{\mathrm{peak}}$ is the peak amplitude of the detected pulse. These are then used to estimate the flux density scales at multiple frequency bands of MWA observations.

As described in Wayth et al. (2017), the performance of the EDA is well characterized, particularly in terms of the receiver and system temperatures, the beam and sky models, and the sensitivity. Further, as shown in their Figure 9, the predicted sensitivity based on the beam and sky models is found to be in good agreement with measurements. Therefore, it is justifiable to use the EDA beam model to calculate the effective collecting area A_e and the system temperature ${T}_{\mathrm{sys}}$, and consequently, the sensitivity A_e/${T}_{\mathrm{sys}}$, or $\mathrm{SEFD}=2{k}_{{\rm{B}}}{T}_{\mathrm{sys}}/{A}_{e}$, for the specific pointing directions of PSR J0437−4715 and PSR J2145−0750. The system temperature ${T}_{\mathrm{sys}}\cong {\eta }_{\mathrm{rad}}{T}_{\mathrm{ant}}+{T}_{\mathrm{rcv}}$, where ${T}_{\mathrm{rcv}}$ is the receiver temperature and ${T}_{\mathrm{ant}}$ is the antenna temperature, calculated as the beam-weighted average sky temperature, using the well-known sky model of Haslam et al. (1982) at 408 MHz (see Equation (3) in Wayth et al. 2017), scaled to the observing frequency as $T{(f)={T}_{408}({f}_{\mathrm{MHz}}/408)}^{\beta }$, where β = −2.55, f_MHz is the frequency (in MHz), and T₄₀₈ is the sky temperature at 408 MHz. As all our EDA observations were made near the meridian transit, this provides fairly reliable estimates of the sensitivity. The values of G and ${T}_{\mathrm{sys}}$ obtained in this manner are then used to estimate the flux density scales at multiple frequency bands of EDA observations (i.e., in 5.12 MHz sub-bands, across the EDA band), in a manner similar to that for MWA pulsar detections (see Section 4.2.1).¹⁶

MSPs generally show complex pulse profiles, with a large number of components, unlike most normal pulsars (e.g., Dai et al. 2015). The degree of complexity often extends to polarization properties, including highly complex polarization position angle variations across the pulse phase, which cannot be well described by the rotating vector model (e.g., Xilouris et al. 1998). Moreover, the emission is often seen to extend over a significant fraction of the pulsar period (W₁₀ ≳ 0.5P, where W₁₀ is the width of the pulse measured at 10% of the peak flux, and P is the rotation period). Furthermore, the conventional radius-to-frequency mapping (Komesaroff et al. 1970; Cordes 1978), which often manifests as a considerable narrowing of the pulse with an increase in the observing frequency, is seldom seen with MSPs (e.g., Kramer et al. 1998). Overall, MSP profiles have been quite challenging to model in terms of emission geometries, which may be due to multiple emission regions within the magnetosphere and pulse components originating at different locations (see Gil & Krawczyk 1997; Kramer et al. 1998; Xilouris et al. 1998; Dyks et al. 2010).

The recent work of Dai et al. (2015) has been very instructive in understanding the frequency evolution of MSP profiles. Using 6 yr of PTA data from Parkes, they were able to perform phase-resolved spectral studies of MSP emission properties, including polarization variations across the pulse phase. PSR J0437−4715 is known to exhibit rather complex, and quite remarkable, evolution in its mean pulse shape (e.g., Navarro et al. 1997; Bhat et al. 2014). Further, as noted by Dai et al. (2015), the profile features in their high-quality Parkes data are generally consistent with our published MWA observations (at a frequency of ∼200 MHz), where a central bright component is flanked by multiple outer components (see Figure 1). Parkes data at very high signal-to-noise ratios (S/Ns; see Figure A1 of Dai et al. 2015) also reveal the emission extending out to more than 85% of the pulse period (i.e., ∼300° in longitude), which was first reported by Yan et al. (2011; though also seen by Navarro et al. 1997). This is now also seen in our very high quality MWA detection (S/N ∼ 3400) shown in Figure 1 (the right panel), where the detectable emission extends to approximately 75% of the pulse period (i.e., ∼260° in longitude).

A comparison of MWA and Parkes data shows that the profile evolution of this pulsar is largely due to its highly complex spectral index (α) variation across the pulse. As seen from the work of Dai et al. (2015), α varies from −0.8 to −2.2 across the large emission window; while the leftward outer-edge component has a relatively flatter spectrum (α > −1.5), the emission within a ∼0.2 phase range left of the central component shows a relatively steeper spectrum (−1.7 > α > −2.2). The emission within a ∼0.2 phase range right of the central component shows a more complex variation, whereby α varies from −1 to −2.2 rather smoothly, and is marked by two prominent minima and maxima. As seen from Figure A1 in Dai et al. (2015), the leading and trailing outer components have generally steeper spectra, which may explain their rapid evolution at the MWA's frequencies.

The pulsar profile is also known for its unusual double-"notch" feature (see Figure 1; near phase ≈0.62 in the 438 MHz profile), besides a notch-like feature near the center, particularly in observations at frequencies in the 0.4–1.5 GHz range (Navarro et al. 1997). It appears that the disappearance of the central notch-like feature may be due to the complex spectral index changes. For instance, the rather abrupt change of α from −2 to −1 closer to the center (i.e., near phase 0 in Dai et al. 2015) may explain the ultimate low-frequency dominance of one component and disappearance of this central notch-like feature.

The disappearance of the double-notch feature in MWA data is likely due to the limitation in our temporal resolution. The width of this feature is approximately 0.01 phase (Navarro et al. 1997), i.e., ∼50 μs, and hence as such is difficult to resolve with our 100 μs native time resolution.¹⁹ The location of this notch feature also happens to be near a local minimum in the spectral index variation, and the resulting nonuniform evolution of the components could also potentially lead to its disappearance in the MWA band. On the other hand, the leftward outer-edge region of the profile (i.e., phase <0.3 in Figure 1) exhibits a relatively flatter spectrum (α > −1.5), which may also explain its reduced prominence at the MWA's frequencies.

PSR J2145−0750 presents yet another interesting case. As evident from Figure 3, the relative strengths of its leading and trailing components continue to decrease at lower observing frequencies and are reversed in the MWA band. Such a reversal was noted earlier by Kuzmin & Losovsky (1996) in their observations at 102 MHz, and more recently in LOFAR data reported in Kondratiev et al. (2016). This can be readily explained given the results of Dai et al. (2015), whose data we re-present in Figure 3. As can be seen from their Figure A23, the leading component of the pulse has α within the range −1.5 > α > −2, whereas for the trailing component it is −2 > α > −2.5. This readily explains the observed reversal of the relative strength and its rapid evolution within the MWA band (Figure 3, left panel). The high time resolution of Parkes data also reveals finer structure, which is not resolved in MWA observations, largely due to the limited time resolution of VCS (100 μs) and non-negligible dispersive smearing (from our 10 kHz channels). Further, the precursor component that appears at 0.2 phase offset prior to the leading peak is not seen in MWA data; this is quite perplexing considering its expected amplitude of ∼10%–40% of the leading peak, assuming the relatively steep spectral index (α ∼ −2) that was inferred by Dai et al. (2015). This may suggest either a possible turnover of the precursor emission or its plausible intermittent nature.

The dynamic spectrum $S(f,t)$ is a 2D record of pulse intensity in time and frequency. It is the most basic observable for scintillation analysis. Figures 6 and 7 show such spectra for PSR J0437−4715 and PSR J2145−0750 obtained over 0.5−1 hr durations, over a continuous 30.72 MHz bandwidth (Figure 6), as well as 9 × 2.56 MHz bands for PSR J2145−0750 (Figure 7). Diffractive scintillation is seen as rapid, deep modulation of pulse intensity in time and frequency, whereas the drifting of intensity maxima (in the time–frequency plane) is generally attributed to refraction through the ISM. Drifting scintles can also be explained by the relative velocities of Earth and the pulsar projected onto the scattering screen; interference with the Doppler-shifted scattered radio waves causes periodic temporal variability of pulse intensity. These temporal variations give rise to drifting when combined with periodic spectral variability due to the temporal delay of the scattered rays (Walker et al. 2004). Drifting scintles are most readily seen in observations of PSR J0437−4715 (Figure 6), where the dynamic spectrum is dominated by a single bright scintle that extends over a timescale longer than our observing duration. As seen from these figures, scintles vary in their brightness, size (in both time and frequency), and orientation (i.e., drift rates); however, average properties can be meaningfully characterized by computing a 2D ACF of $S(f,t)$ in frequency and time lags: $\rho (\nu ,\tau )=\langle S(f,t)$ $\,S(f+\nu ,t\,+\tau )\rangle$, where ν and τ are the frequency and time lags, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The ACFs computed in this manner are shown in the bottom panels of Figure 6. The characteristic widths in frequency and time, i.e., scintillation bandwidth ${\nu }_{{\rm{d}}}$ and timescale ${\tau }_{\mathrm{iss}}$, respectively, can be determined by fitting a 2D elliptical Gaussian function to these ACFs. Following the published literature (Gupta et al. 1994; Bhat et al. 1999; Wang et al. 2005), we fit a functional form ${\rho }_{g}\,(\nu ,\,\tau )={C}_{0}$ $\,\exp [-({C}_{1}{\nu }^{2}+{C}_{2}\nu \tau \,+{C}_{3}{\tau }^{2})]$, to yield characteristic scales: the decorrelation bandwidth, ${\nu }_{{\rm{d}}}={(\mathrm{ln}2/{C}_{1})}^{0.5}$, measured as the half-width at half-maximum of the correlation peak; scintillation timescale, ${\tau }_{\mathrm{iss}}={(1/{C}_{3})}^{0.5}$, measured as the 1/e width of the correlation maximum; and the drift rate, ${dt}/d\nu =-({C}_{2}/2\,{C}_{3})$. In addition, a "drift-corrected" scintillation bandwidth can also be estimated, i.e., ${\nu }_{{{\rm{d}}}_{{\rm{c}}}}\,={(\mathrm{ln}2)}^{0.5}\,{({C}_{1}-{C}_{2}^{2}/4{C}_{3})}^{-0.5}$, to reduce the underestimation bias resulting from significant refraction typically seen in our data. The parameters determined in this manner are tabulated in Table 3, along with some physical properties that can be derived, namely, the mean turbulence strength $\overline{{C}_{{\rm{n}}}^{2}}$, the strength of scattering u, and the scintillation pattern speed V_iss. Further, the diffractive and refractive scattering angles (${\theta }_{\mathrm{diff}}$ and ${\theta }_{\mathrm{ref}}$) can be estimated from the measurements of decorrelation bandwidth and the drift rate dt/dν, respectively, under the assumption of the thin screen scattering model, which is generally employed in the context of scintillation (see Cordes et al. 1986; Gupta et al. 1994; Bhat et al. 1999).

Table 3.  Scintillation Parameters from MWA Observations

PSR	${\nu }_{{\rm{d}}}$	${\tau }_{\mathrm{iss}}$	dt/dν	Viss	${V}_{\mu }$	$\overline{{C}_{{\rm{n}}}^{2}}$	${\theta }_{\mathrm{diff}}$	$| {\theta }_{\mathrm{ref}}|$	sda
	(MHz)	(s)	(${\rm{s}}\,{\mathrm{MHz}}^{-1}$)	(km s−1)	(km s−1)	(${{\rm{m}}}^{-20/3}$)	(mas)	(mas)
J0437−4715	1.43 ± 0.4	613 ± 165	−181	114 ± 36	100	8 × 10−5	0.8	0.2	0.58
J2145−0750	0.027 ± 0.004	300 ± 45	1080	99 ± 18	33	1.4 × 10−4	3	0.2	0.09

Note.

^aThe sparseness parameter that quantifies the time–frequency occupancy of scintles (see the text for details).

Download table as:  ASCIITypeset image

For PSR J0437−4715, we measure ${\nu }_{{\rm{d}}}=1.4\pm 0.4$ MHz,²⁰ which is consistent with that reported in an earlier publication (Bhat et al. 2014). As we summarized in that paper, there are striking discrepancies between previously published measurements, which may be attributed to the choice of suboptimal observing parameters in much of the earlier work. The mean turbulence strength, $\overline{{C}_{{\rm{n}}}^{2}}=8\times {10}^{-5}$ ${{\rm{m}}}^{-20/3}$, is the second lowest (after PSR B0950+08; see Phillips & Clegg 1992) measured among all known pulsars, making PSR J0437−4715 one of the most weakly scattered pulsars. Any timing perturbations resulting from scattering can therefore be safely ignored at its timing frequencies (∼1–2 GHz).

The measured drift rate dt/dν is almost twice as large as our previous observations and of opposite sign, though this can be attributed to time-varying refraction that is commonly seen toward many of the low-DM pulsars (Bhat et al. 1999; Wang et al. 2005). The refractive scattering angle ${\theta }_{\mathrm{ref}}$ ∼ 0.2 mas is of comparable magnitude still, and it is the combination of a large ${\nu }_{{\rm{d}}}$, the pulsar's high proper motion (${V}_{\mu }$ = 100 km s⁻¹), and a long timescale (${\tau }_{\mathrm{iss}}$ = 613 ± 165 s) that gives rise to visibly prominent drift patterns that we see in the MWA band. Such prominence of refraction would generally imply the existence of discrete density structures along the line of sight, which can also give rise to parabolic scintillation arcs that we discuss later in Section 5.2.2. Variability of the measured drift rate can also be explained by the orbital motions of Earth and the pulsar; as their relative velocities projected onto the thin screen vary with time, so does the Doppler shift of radio waves scattered off asymmetric structures in the scattering screen.

For the PSR J2145−0750 observations presented in Figure 6 (over a 30.72 MHz bandwidth centered at 154.24 MHz), we measure ${\nu }_{{\rm{d}}}=27\pm 4$ kHz and ${\tau }_{\mathrm{iss}}$ = 300 ± 45 s. Thus, the scintles are barely resolved (in frequency) with our 10 kHz spectral resolution with the VCS. The derived value for $\overline{{C}_{{\rm{n}}}^{2}}$ is nearly twice that of PSR J0437−4715, making the sight line of PSR J2145−0750 comparatively less anomalous.

The scintillation velocity V_iss can be derived from the measurements of ${\nu }_{{\rm{d}}}$ and ${\tau }_{\mathrm{iss}}$, e.g., in the simplest case where the scattering medium is approximated as a thin screen located at a certain distance between the pulsar and the observer, V_iss is given by (see Gupta et al. 1994)

$$\begin{eqnarray}&&{V}_{\mathrm{iss}}={A}_{\mathrm{ISS}}\,{({\nu }_{{\rm{d}}}Dx)}^{1/2}{({\nu }_{\mathrm{obs}}{\tau }_{\mathrm{iss}})}^{-1},\end{eqnarray} \tag{ 1 }$$

where D = ${D}_{\mathrm{os}}$ + ${D}_{\mathrm{ps}}$ is the pulsar distance, and $x={D}_{\mathrm{os}}/{D}_{\mathrm{ps}}$, i.e., the ratio of the distances from the screen to the observer (${D}_{\mathrm{os}}$) and from the screen to the pulsar (${D}_{\mathrm{ps}}$). The constant, A_ISS, relates the timescale ${\tau }_{\mathrm{iss}}$ to the velocity, for which Cordes & Rickett (1998) derive a value of 2.53 × 10⁴ km s⁻¹ for a Kolmogorov turbulence spectrum and homogeneously distributed medium (with ν_d in MHz, ν_obs in GHz, τ_iss in s, and D in kpc), and Gupta et al. (1994) derive A_ISS = 3.85 × 10⁴ km s⁻¹ for a single asymmetrically located thin screen (and is the value that we adopt for our estimation).

The estimated scintillation velocity V_iss = 99 ± 18 km s⁻¹ (including the uncertainty in the pulsar's distance, D = 530 ± 60 pc) is three times larger than the measured proper motion of 33 km s⁻¹ for this pulsar (Reardon et al. 2016), which may suggest an asymmetric location for the underlying scattering screen.

The 2D power spectrum of the dynamic spectrum S(f, t) is referred to as the secondary spectrum, ${S}_{2}({f}_{\nu },{f}_{{\rm{t}}})\,=\,| {S}^{\dagger }(\,f,t){| }^{2}$ (where † indicates 2D Fourier transform), and is a powerful technique that captures interference patterns produced by different points in the image plane (e.g., Stinebring et al. 2001; Walker et al. 2004; Cordes et al. 2006). The "fringe rates" (in time and frequency), ${f}_{{\rm{t}}}$ and ${f}_{\nu }$, are essentially the fringe frequency and delay parameters, where ${f}_{\nu }$ is a measure of the differential time delay between pairs of rays, and ${f}_{{\rm{t}}}$ is the temporal fringe frequency. Interference between scattered wave fronts from the origin and pairs of points along an axis in the direction of the net velocity vector ${V}_{\mathrm{eff}}$ produces parabolic scintillation arcs, which can be represented by ${f}_{\nu }=\eta \,{f}_{{\rm{t}}}^{2}$. The parabolic arcs are thus essentially a natural consequence of small-angle forward scattering. The fringe frequency and delay parameters ${f}_{{\rm{t}}}$ and ${f}_{\nu }$ can be related to the curvature of the arc (η); following Cordes et al. (2006),

$$\begin{eqnarray}&&\eta =\displaystyle \frac{D\,s\,(1-s)\,{\lambda }^{2}}{2\,c\,{V}_{\mathrm{eff}}^{2}\,{\cos }^{2}\theta },\end{eqnarray} \tag{ 2 }$$

where D is the pulsar distance, s is the fractional distance of the screen from the source (i.e., $s={D}_{\mathrm{ps}}/D$, and s = 0 at the source), ${D}_{\mathrm{eff}}\equiv D\,s\,(1-s)$ is the effective distance to the screen, and θ is the angle between the net velocity vector V_eff and the orientation of the scattered image. The quantity ${V}_{\mathrm{eff}}$ is the velocity of the point in the screen intersected by a straight line from the pulsar to the observer, and it is the weighted sum of the pulsar's binary and proper motions (${V}_{\mathrm{bin}}$ and ${V}_{\mu }$, respectively) and the motion of the screen and the observer (${V}_{\mathrm{scr}}$ and ${V}_{\mathrm{earth}}$, respectively). Its transverse component V_eff⊥ (see Equation (3) in Bhat et al. 2016) determines the measured timescale ${\tau }_{\mathrm{iss}}$. Thus, the measurement of η can be used to determine the location of the scatterer, when all the contributing terms are precisely known.²¹

In a previous publication (Bhat et al. 2016) that reported our first observations of scintillation arcs in PSR J0437−4715, we used the measurement of η and the knowledge of the pulsar's precise orbital and astrometric parameters, to determine the location of the underlying scatterer to be at 115 ± 3 pc. Interestingly, this was consistent with the location inferred from independent observations from Parkes made 2 weeks later at a frequency of 732 MHz.

Figure 8 shows the secondary spectrum from our new observations of PSR J0437−4715 (for the dynamic spectrum shown in Figure 6); the resolutions in fringe frequency and delay axes are 0.195 mHz and 0.195 μs, respectively (i.e., identical to those in Bhat et al. 2016). The parabolic arc features are consistently seen across different frequency segments of these new observations, albeit somewhat less prominently in the upper one-sixth of the band, where the pulsar is visibly dimmer compared to the remainder of the observing band. Furthermore, the arc feature is also visibly stronger than our previous detection, albeit more diffuse in appearance, which is reminiscent of a filled parabola (e.g., Putney & Stinebring 2006). Following a feature-extraction technique as described in Bhat et al. (2016), we compute an "arc strength" parameter ${P}_{\mathrm{arc}}$, given by

$$\begin{eqnarray}&&{P}_{\mathrm{arc}}(\eta )=\displaystyle \frac{1}{N}\sum _{i=1}^{N}{S}_{2}(\eta {f}_{{\rm{t}},{\rm{i}}}^{2},{f}_{{\rm{t}},{\rm{i}}}).\end{eqnarray} \tag{ 3 }$$

This method, which can be compared to the 1D generalized Hough transform, is expected to be more robust, particularly when the arc feature is limited by S/N, which is the case with our detections. The summing procedure is performed along the points of arc outside an excluded low-frequency noise region, and out to delays beyond which little power is detectable (i.e., i = 1...N, where N corresponds to ${f}_{\nu }$ = 20 μs). Further, in order to better account for the visibly larger diffused arc feature, we consider a slightly modified parameterization for the "thickness" of the parabola, ${p}_{{\rm{d}}}$, which is essentially the distance (in pixel units) from the arc curvature that is used for the computation of S₂, e.g., ${p}_{{\rm{d}}}=2$ means that only those spectral points lying within 2 pixels of the parabola (either horizontal or vertical) are included in the sum. This is expected to yield a more robust parameterization, particularly when the arc feature is more diffuse in nature.

Zoom In Zoom Out Reset image size

The right panel of Figure 8 shows our estimated P_arc(η) for the secondary spectrum shown in the left panel of this figure. The curve appears to be more broadly peaked than that of our earlier detection, despite its visibly higher degree of prominence. A closer examination of the secondary spectrum also reveals that the arc feature is more prominent in the right half. The best-fit value of the curvature of the arc η, which corresponds to the maximum value of ${P}_{\mathrm{arc}}$, is η = 0.79 ± 0.03, for the case where ${p}_{{\rm{d}}}=1.5$; however, the estimate is somewhat more uncertain if we consider the more broadly peaked curve for ${p}_{{\rm{d}}}=2.5$ (i.e., a thicker parabola). Nevertheless, all three curves peak at similar values of η.²² The implied distance to the screen (after accounting for the expected ${V}_{\mathrm{eff}}$ at our observing epoch; for details, see Bhat et al. 2016) is s = 0.26 ± 0.02 (or 115 ± 3 pc), which is in excellent agreement with our previously reported measurements. This stability of the η value reaffirms the presence of a compact scatterer in the line of sight to the pulsar.

As seen from Equation (2), the inferred location s depends also on the orientation θ of the scattered image relative to V_eff (and is generally an unknown). However, since the 3D sky geometry of the pulsar is precisely known (van Straten et al. 2001), it is indeed possible to calculate the projection of the vector velocity V_eff relative to the scattering screen (assuming that it is on the plane of the sky). The effective velocity is still dominated by the proper motion of the binary system, and for the two epochs of MWA observations, incidentally, the projection angles differ by only 26. Therefore, any resultant changes in η are not easily measurable given our uncertainties.

Simultaneous multifrequency observations such as those presented in this paper can also be used, in principle, to investigate the frequency scaling of scintillation parameters. Theoretical treatments based on a Kolomogorov-type turbulence spectrum for electron density fluctuations in the ISM predict ${\nu }_{{\rm{d}}}\propto {\nu }^{4.4}$ and ${\tau }_{\mathrm{iss}}\propto {\nu }^{1.2}$ (e.g., Cordes et al. 1986). However, departures from these theoretical scalings have been seen toward a number of objects, and on a global scale, observations indicate a shallower scaling (e.g., ${\nu }_{{\rm{d}}}\propto {\nu }^{3.9\pm 0.2};$ Bhat et al. 2004). An important caveat is that, at low frequencies, long-term (refractive) modulations in ${\nu }_{{\rm{d}}}$ and ${\tau }_{\mathrm{iss}}$ can bias the measured value at a given epoch of observation (e.g., Gupta et al. 1994; Bhat et al. 1999), and ideally observations spanning many refractive cycles are needed for a reliable estimation of the frequency scaling index.

Notwithstanding this, we attempted the related analysis using our multiband observations. PSR J0437−4715 data were not suitable for this purpose, as the expected ${\nu }_{{\rm{d}}}\gtrsim 1.28$ MHz bandwidth for a large subset of our data (at frequencies ≳150 MHz). However, our measurement of ${\nu }_{{\rm{d}}}\sim 35$ kHz for PSR J2145−0750 in the 140–170 MHz band (see Figure 6) suggests that ${\nu }_{{\rm{d}}}$ is measurable within the constraints of our 10 kHz spectral resolution (and 2.56 MHz sub-bands) at frequencies ≳140 MHz.

Figure 7 shows the dynamic spectra from our PSR J2145−0750 observations shown in Figure 2. As seen from this figure, scintles are resolved within our observing parameters, i.e., 10 kHz ≲ ${\nu }_{{\rm{d}}}$ ≲ 2.56 MHz over the 112–220 MHz frequency range of our observation. This demonstrates the value of the MWA's observable frequency range for studying scintillation parameters as a function of frequency. However, because of the limited S/N, a full 2D correlation analysis proved less meaningful, and we therefore adopted a relatively simpler analysis, whereby we make use of 1D cuts of the 2D ACF along ν = 0 and τ = 0 axes to estimate ${\nu }_{{\rm{d}}}$ and ${\tau }_{\mathrm{iss}}$ (by fitting for a Gaussian function).

Our estimates obtained in this manner are shown in Figure 9, and the derived scaling indices are 2.3 ± 0.3 for ${\nu }_{{\rm{d}}}$ and 0.5 ± 0.3 for ${\tau }_{\mathrm{iss}}$, i.e., much shallower than we expect based on the empirically determined scaling or the theoretical predictions. However, this is perhaps less surprising, as our estimates are derived from a single epoch of observations and the refractive modulation may not necessarily be correlated across our large fractional bandwidth. This is probably best alleviated by multiple similar observations that span a large number of refractive cycles and by employing more optimal observing strategies (e.g., fewer sub-bands with larger sub-bandwidth to allow 2D correlation analysis and more reliable estimates).

Zoom In Zoom Out Reset image size

The dynamic spectra of PSR J2145−0750 in Figure 7 are marked by a visible sparseness in the distribution of scintles. This time–frequency occupancy of bright features is likely related to the filling fraction, f_d, typically used in the context of scintillation in order to account for the density (or a finite number) of scintles in the dynamic spectrum. In the strong scattering regime, bright features are expected as a natural consequence of an exponential distribution of intensity variations, with separations in time and frequency much larger than the characteristic sizes of scintles (Cordes 1986). Observational estimations of scintillation parameters have adopted values in the range 0.2 < f_d < 0.5 (Gupta et al. 1994; Bhat et al. 1999) and so are larger than a more conservative value of 0.01 assumed by Cordes (1986). Observations show that this filling fraction tends to vary between pulsars. In order to quantify this time–frequency occupancy of scintles in our data, we introduce the following "sparseness" parameter:

$$\begin{eqnarray}&&{s}_{{\rm{d}}}=\displaystyle \frac{1}{{N}_{{\rm{f}}}\,{N}_{{\rm{t}}}}\sum _{i=1}^{{N}_{f}}\sum _{j=1}^{{N}_{t}}{s}_{\mathrm{on}}({f}_{i},{t}_{j}),\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray*}{s}_{\mathrm{on}}({f}_{{\rm{i}}},{t}_{{\rm{j}}})=\left\{\begin{array}{ll}1 & \mathrm{if}\,S({f}_{{\rm{i}}},{t}_{{\rm{j}}})\gt 5{\sigma }_{\mathrm{off}}({f}_{{\rm{i}}},{t}_{{\rm{j}}})\\ 0 & \mathrm{otherwise}\end{array}\right.,\end{eqnarray*}$$

where σ_off(f_i, t_j) is the off-pulse rms corresponding to the dynamic spectral point S(f_i, t_j). In other words, this represents the conservative case, where at least one sample within the on-pulse window contributing to S(f_i, t_j) exceeds 5σ above the mean off-pulse level (and hence some detectable pulse energy). In some ways, this can be treated as an upper limit to the filling fraction.

For PSR J2145−0750 data shown in Figure 7, values of the s_d parameter computed in this manner range from ∼0.05 to ∼0.12, the lower values being measured for the highest two frequency bands. It is possible that the relatively low S/Ns of these data may be somewhat biasing the estimation of s_d. However, for the data shown in Figure 6, which are limited neither by S/N nor by small-number statistics in terms of the number of scintles, the above estimation still yields s_d ∼ 0.09.²³ This is comparable to the values for the data shown in Figure 7 and is significantly lower than that estimated for the dynamic spectra of PSR J0437−4715 in Figure 6, for which we estimate s_d ∼ 0.58. It is possible that such low levels of time–frequency occupancy in terms of scintles are a characteristic of certain lines of sight, where underlying small-scale density structures deviate from the standard Kolmogorov-type distribution. A more quantitative analysis is deferred to a future paper as further observations accrue from our ongoing monitoring project.

Using the procedures detailed in Section 4.2, we estimated the flux density measurements for the data presented in this paper. Since our MWA data are limited in both S/N and temporal resolution (besides multiple caveats of insufficient averaging of flux modulations due to scintillation effects), we limit ourselves to a fairly simple analysis, where we ignore the spectral index variation across the pulse profile, and consider just the spectrum of the average emission (i.e., integrated pulse profiles as a whole). Figure 10 presents our measured flux densities for PSR J0437−4715 and PSR J2145−0750, along with the published values from Dai et al. (2015) and Dowell et al. (2013). Parkes measurements are from regularly sampled Parkes PTA data over a 6 yr time span and therefore are more reliable; however, we note that both pulsars are seen to be highly variable even in the Parkes bands, e.g., rms variability ∼50%–70% is noted for PSR J0437−4715, whereas PSR J2145−0750 is seen to vary much more, with rms variation comparable to mean flux densities in the 50 and 10 cm data (i.e., 732 MHz and 1.4 GHz). Similarly, large variations are also seen in our MWA observations, and therefore the implied spectral behavior can be considered only indicative. EDA measurements are more reliable than those from the MWA, as they are obtained from multiple observations spanning several months; for PSR J2145−0750, nine observations over a 4-month time span were used, whereas six observations spanning ∼8 months were used in the case of PSR J0437−4715. Still, this only sparsely samples long-term refractive cycles, and an insufficient averaging of resultant modulations in flux densities can potentially skew our estimates. Table 4 presents a summary of our spectral index estimates.

Zoom In Zoom Out Reset image size

For PSR J2145−0750, excluding LWA measurements, the combination of EDA and Parkes measurements yields a spectral index estimate, α_e = −1.73 ± 0.07; however, a larger value of α_m = −2.12 ± 0.07 is estimated if we use MWA measurements instead. The reported value of α_p = −1.96 ± 0.06 by Dai et al. (2015) is still within this range and agrees with our estimates at the 2σ level, which is encouraging and reaffirms a generally steeper spectrum of this pulsar. Furthermore, it appears that our measurements (in particular, those from the EDA) are not consistent with a suggested spectral turnover for this pulsar by Dowell et al. (2013). In fact, their quoted flux densities are based on a limited number of observations and hence are subject to possible bias from insufficient averaging of long-term refractive modulations, a caveat also applicable for MWA and EDA measurements. Therefore, the possibility of temporal variations (in flux densities) causing the observed inconsistency cannot be precluded.

For PSR J0437−4715, a number of factors prevented us from obtaining reliable estimates of its flux densities. It was difficult to accurately calibrate MWA observations shown in Figure 1 for reliable flux scales because a large number of tiles (58 out of 128) were flagged during the calibration procedure. The measurements shown in Figure 10 are hence from observations shown in Figure 6(a) (i.e., over 140–170 MHz, broken up into 12 × 2.56 MHz contiguous bands). Furthermore, this pulsar also shows substantial variability in apparent brightness—as much as by a factor of ∼6 between different observations. Given these, the indicated α_e = −1.1 ± 0.07 (from EDA and Parkes data), which is obtained by excluding the observing epochs when the pulsar was significantly brighter than average, can only be treated as an indicative upper limit (i.e., $| \alpha | \gt 1.1$). The spectral estimate from MWA and Parkes data is somewhat steeper (α_m = −1.41 ± 0.14); however, it is still shallower than the value (α_p = −1.67 ± 0.04) from Dai et al. (2015) from Parkes-only measurements. We further note that the earlier work of McConnell et al. (1996), which reported the first detection of this pulsar at frequencies below 100 MHz, was inconclusive in terms of a possible turnover at low frequencies. Further observations are needed to explore this in detail.