MeerKAT Observations of the Reversing Drifting Subpulses in PSR J1750-3503

The data used in this paper were recorded with the MeerKAT telescope using the PTUSE ⁵ single-pulse observing mode described in Bailes et al. (2020) and were collected during seven observing sessions (see Table 1) with the L-band receiver. The first two observing sessions were taken using all available MeerKAT antennas, while the following sessions were taken in subarray mode to optimize the telescope use (see Song et al. 2021). The central frequency of the observations is 1284 MHz, while the bandwidth is 642 and 856 MHz for data collected in 2019 and 2020, respectively, and at a 38.3 μs time resolution (see Table 1 for more details). The number of frequency channels were 768 and 1024 for the observations in 2019 and 2020, respectively, corresponding to a common channel width of 0.836 MHz. The data were processed using the MeerTime Single Pulse pipeline (Keith et al., in preparation). The pipeline produces science-ready data using three stages of processing. In the first stage, the time-series data from the telescope are folded into a single PSRFITS archive using dspsr (van Straten & Bailes 2011). Polarization calibration is applied, and the frequency channels contaminated by radio frequency interference are identified using noise statistics. In the second stage of processing, single-pulse data are formed, after applying the polarization calibration and channel mask obtained from the first stage. The pipeline outputs one set of single-pulse archives that are averaged over the full bandwidth and one set that is frequency scrunched to 16 subbands. The third stage of the pipeline produces a number of diagnostic plots to assess the data quality. More details of this pipeline will be presented elsewhere (Keith et al., in preparation). For the work presented here, we have used the band-averaged single-pulse archives produced by this pipeline.

In Figure 1 we show sequences of single pulses for all observing sessions. Except for the first session (the shortest), the collected single pulses are characterized by a very unstable drifting behavior (see panels (b)–(g)).

Zoom In Zoom Out Reset image size

In Figure 2 we show the pulse energy distributions for 1031 pulses of the second observing session as two histograms corresponding to the on and off pulse energies. The off-pulse distribution is centered around zero and reflects the noise characteristics of the baseline level. The lack of a bimodal shape in the on-pulse curve indicates that there is no obvious evidence for nulling.

Zoom In Zoom Out Reset image size

Constraining the viewing geometry of PSR J1750−3503 is key to applying the carousel model and understanding how it gives rise to the observed pattern of drifting subpulses. This geometry is parameterized by two quantities: the magnetic inclination angle α and the impact angle β of the observer's line of sight. A commonly used method to determine these is studying the pulsar polarization.

Figure 3 shows the integrated profile of PSR J1750−3503 and its polarization properties. In the top panel, total intensity is shown by the black curve, the linearly polarized fraction by the red curve, and the circularly polarized fraction by the green curve. The lower panel shows how the position angle (P.A.) of the linear polarization varies as a function of pulse longitude.

Zoom In Zoom Out Reset image size

There is no evidence of a correlation between the drifting subpulses and polarization properties such as is seen in some pulsars when analyzing the individual pulses (e.g., PSR B0031−07; Ilie et al. 2020).

The P.A. curve was analyzed with PSRSALSA ⁶ (Weltevrede 2016) and is shown in Figure 3. It is perfectly compatible with the rotating vector model (Radhakrishnan & Cooke 1969), although fitting this was not sufficient to constrain α. However, the steepest gradient of the curve was −2.2 ± 0.3, which is equal to sinα/sinβ (Komesaroff 1970). Better constraints can be placed on α and β by taking into account the large observed profile width W. The width of the profile at 10% of its peak intensity is W₁₀ = 58° ± 7°, found by fitting a model composed of two von Mises functions. The pulse profile width that is observed depends on how the line of sight intercepts the conical emission beam, as determined by the viewing geometry α and β. Gil et al. (1984) showed that the profile width W can be determined by

$\begin{eqnarray}&&\cos \rho =\cos \alpha \cos (\alpha +\beta )+\sin \alpha \sin (\alpha +\beta )\cos \left(\displaystyle \frac{W}{2}\right),\end{eqnarray} \tag{ 1 }$

where ρ is the half-opening angle of the cone of radio emission. Radio emission of a given frequency is assumed to be produced at some altitude h_em in the magnetosphere, within the region delimited by the tangents to the last open field lines touching the light cylinder, which form a conical beam. The half-opening angle ρ is given by

$\begin{eqnarray}&&\rho \approx \sqrt{\displaystyle \frac{9\pi {h}_{\mathrm{em}}}{2{cP}}}\end{eqnarray} \tag{ 2 }$

in the small-angle limit (h_em ≪ r_LC, e.g., Rankin 1990), where c is the speed of light. The emission height of radio pulsars is believed to lie within 200 to 400 km (Johnston & Karastergiou 2019) irrespective of period. This means that for PSR J1750−3503, where P = 0.684 s, the cone opening angle lies in the range 5° ≲ ρ ≲ 10° if the emission height is not atypical.

Combining Equations (1) and (2) leads to a set of constraints such that only certain combinations of α and β are "allowed." These allowed geometries indicate that ∣β∣ < 10° and 5° ≲ α ≲ 25°. These values are input for the subpulse-drift interpretation and modeling that we pursue next.

To study the drift characteristics in J1750−3503, we use longitude-resolved fluctuation spectra (LRFS, Backer 1970), whose computation involves discrete Fourier transforms of consecutive pulses along each rotational phase longitude. To find the time intervals in which the drifting behavior is stable, we calculate the LRFS of intervals of 128 single pulses, shifting the starting point by one period every time (see also Serylak et al. 2009). The repetition time of the drift pattern, P₃, is determined by fitting a Gaussian near the major peak in the longitude-averaged fluctuation spectra. In Figure 4 we show the variation of the fluctuation spectra as a function of the pulse number corresponding to the interval start. Visual inspection of the drift in Figure 1 shows the patterns are not always stable but vary. The continuous range of these periodicities is visible in the bottom panels of Figure 4. Below we focus on a number of notable periodicity epochs within each session. The repetition time seems to be stable only for the first observing session, with P₃ = (43 ± 1)P. To calculate the repetition time for the first session with a better accuracy we use the PSRSALSA software package (see Weltevrede (2016), for more details), which results in P₃ = (43.5 ± 0.4)P. The measurements for the other observing sessions are characterized by much more variability, resulting in session-averaged repetition times P₃ = (47 ± 14)P, P₃ = (43 ± 8)P, P₃ = (44 ± 6)P, P₃ = (47 ± 7)P, P₃ = (50 ±13)P, and P₃ = (40 ± 3)P, respectively, for the observing sessions presented in panels (b)–(g). The errors on these mean values are the 1σ standard deviation. They signify the variation in the underlying pattern. During those sessions, a few intervals with a stable P₃ can be identified (see Figure 4). These show steady repetition values of P₃ = (43 ± 1)P, P₃ = (41 ± 1)P, P₃ = (35 ± 1)P, P₃ = (41 ± 1)P, P₃ = (37.8 ± 0.3)P, P₃ =(41 ± 1)P, and P₃ = (39 ± 2)P, respectively, for the starting pulse numbers 1–167 in panel (a), 660–773 and 852–903 in panel (b), 1–88 in panel (e), 8–73 in panel (f), and 20–103 and 228–318 in panel (g). Furthermore, a number of monotonically increasing and decreasing P₃ values are identified. In panel (b), between pulse numbers 131 and 165 the repetition time P₃ decreases from (61 ± 3)P to (44 ± 1)P, and between 463 and 534, P₃ increases from (39 ± 2) to (46 ± 2)P. Between the starting pulse 64–158 in panel (d) P₃ increases from (34 ± 1)P to (46 ± 4)P. In panel (f), between pulses 74 and 115, P₃ increases from (41 ± 2)P to (52 ± 1)P, and between pulses 135 and 156, it increases from (49 ± 2)P to (72 ± 5)P.

Zoom In Zoom Out Reset image size

We identify seven sections with relatively stable P₃ and five sections with monotonically increasing/decreasing P₃, while no evidence of periodic changes in the measured P₃ values was found. Only during the first observing session is the duration when P₃ is stable larger than the length of LRFS. Except for the first observing session, inspection of single pulses (see Figure 1) shows that time intervals with both stable and changing P₃ are characterized by drift direction changes. Since the analysis performed using the LRFS is limited by the length of the Fourier transformation, we perform the drift analysis using single pulses in the following section.

To find the subpulse location (see the red crosses in the upper panels of Figure 5), the data of a single pulse are convolved with a Gaussian function having the mean width of the subpulses. The peaks in the convolution are used to determine the longitudes of the subpulses. In the analysis presented in this paper, we consider only subpulses with a signal-to-noise ratio higher than 5. Furthermore, the detected peaks are visually inspected and only subpulses that are located together as a drift band, and are within the expected longitude range, are considered in further analysis. In Figure 6 we show the detection procedure for a sample single pulse, which results in the following detection fractions for all observing sessions: 70%, 83%, 69%, 74%, 69%, 76%, 79%, that is, data presented in panels (a)–(g), respectively, in Figures 1 and 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.1. Drift Rate

3.2.2. Subpulses Separation P₂

Another quantity that we can use to characterize the drifting subpulses is the longitudinal separation between subpulses, P₂, or in other words, the separation between adjacent drift bands. We use the drift rate information and subpulses locations in single-pulse data to analyze the subpulses separation. Figure 7 shows the distribution of P₂ measurements for all single pulses, P₂ = (18.58 ± 0.07)°, single pulses with positive drift, P₂ = (18.79 ± 0.07)°, and single pulses with negative drift P₂ = (17.52 ± 0.16)°. The uncertainties correspond to the standard error.

Zoom In Zoom Out Reset image size

Previous studies have shown that for a number of pulsars a variation of drift rate is associated with the changing average profile shape (see, e.g., Esamdin et al. 2005, and references therein). To show the variation of the average profile in different drift states, in Figure 8 we compare the average profiles produced by single pulses with positive and negative drift rates. The profiles were normalized using the number of single pulses used. The black solid line in the figure corresponds to the average profile for all single pulses. The average profile for positive drift is characterized by one skewed peak formed from the underlying two-peak structure. On the other hand, in the average profile for single pulses with negative drift, the two-peak structure is more visible. Moreover, the maximum peak seems to be shifted toward earlier longitudes. To estimate differences in intensity, we calculate areas under the curves and normalize them assuming that the area under the average profile for all single pulses is equal to one. The area for positive and negative drift, respectively, is 1.01 ± 0.06 and 0.98 ± 0.06, and here the uncertainties correspond to the standard error. The area measurements are consistent within the uncertainties, suggesting no significant changes in subpulse brightness during reversals.

Zoom In Zoom Out Reset image size

3.3.1. Drift Direction Change

To explore the nature of drift direction change in more detail we use a method based on a linear regression, Segmented. ⁸ Figure 9 shows a detailed analysis of subpulse positions for selected time intervals. We use the coefficient of determination, R², to check how well the data are replicated by the models. For all tracks shown in Figure 9, R² is 0.77 ± 0.11 and 0.75 ± 0.14 for the cubic smoothing spline and linear regression fits, respectively. The colored solid lines in the upper panels of Figure 9 show the linear regression fits, while the bottom panels show a comparison of drift rates derived by the cubic smoothing spline fits (solid blue line) and the linear regression fits (points with error bars). The results suggest that the drift in J1750−3503 can be equally well described by either smooth or rapid changes in drift direction. The vertical dashed lines in the upper panels of Figure 9 mark the estimated time of drift direction change. The break points are consistent within uncertainties throughout all subpulse tracks (see, e.g., pulse number 320 in panel (a), pulse numbers 975 and 1010 in panel (d), and pulse numbers 21 and 50 in panel (h), of Figure 9).

Zoom In Zoom Out Reset image size

The pulsed emission is visible for a short duration every rotation, which could mean that the measured apparent drift rate of drifting subpulses is the result of the aliasing effect. The shortest actual P₃ unaffected by aliasing, that is, the time in which the consecutive subpulse shifts to the space occupied by the previous one, is twice the pulsar rotation period P. The actual modulation period, P₃, can be expressed as P/P₃ = n +P/P₃ ^obs. Here P₃ ^obs is the apparent (possibly aliased) modulation period, where −0.5 < P/P₃ ^obs < 0.5. The sign of P₃ ^obs changes depending on the apparent sign of the gradient of the drift bands. The alias order n is the nearest whole number of undetected subpulses from one pulse to the next (see, e.g., Edwards & Stappers 2002). Note that the definition of n differs from the alias order introduced in Gupta et al. (2004). In order to estimate the true value of P₃ the alias order has to be resolved, but in general it is not always possible to infer the value of n from the observations. In this section we explore the possibility that the peculiar drift characteristics of J1750−3503 are caused by aliasing.

To model single pulses we use the approach proposed in Szary et al. (2020) with solid-body-like rotation of sparks around the center of the polar cap. To achieve the solid-body-like rotation we assumed that the angular velocity of all sparks is equal. At times when the drift is relatively stable, the measured repetition time P₃ ^obs ∼ 44P (see Figure 4 (a)). In Figure 10 we compare single pulses observed during the first observing session (left panel) and the modeled single pulses (right panel). The pulses are modeled using the following pulsar geometry: inclination angle α = 10°, impact parameter β = − 45, and the emission height h_em = 300 km, which results in the opening angle ρ = 83 (see Equations (1) and (2)) and the measured pulse width W₁₀ = 44° ± 2°. To determine the number of subbeams in the carousel, and therefore the sparks at the polar cap, for this specific geometry, we compare phase information obtained using the LRFS of the observed and modeled signals. The single pulses are modeled using P₃ = 0.9778P, resulting in the observed repetition time P₃ ^obs = (43 ± 3)P. We use 13 subbeams with an imprint on the polar cap with the Gaussian variance of 9 m (taking the neutron star radius to be 10 km). Note that since the drifting phenomenon is affected by the viewing geometry (Gupta et al. 2004), the derived number of subbeams is a unique solution only for a given pulsar geometry.

Zoom In Zoom Out Reset image size

At times when the drift in J1750−3503 is characterized by nonstable behavior (see, for instance, panel (b) in Figure 4), P₃ ^obs ranges from ∼33P to ∼85P. Assuming the modulation is observed in the first alias order, and assuming the actual drift direction is the same as that observed (n = 1), the actual repetition time, P₃, ranges from 0.97P to 0.988P. If, on the other hand, we assume that the actual drift direction is opposite to that observed (n = − 1), the actual repetition time, P₃, ranges from 1.011P to 1.031P. As follows from the above calculations, a change in P₃ of only about 6% could produce both the observed drift direction change and the measured P₃ ^obs range. Furthermore, for higher alias orders an even smaller change in P₃ results in the observed drift characteristics. To model the changing P^obs ₃ we alter the angular velocity of sparks according to the measured drift rates (see the bottom panels in Figure 5). In Figure 11 we compare the single pulses observed during the second observing session (top panel) with the modeled single pulses (middle panel). There, P₃ has a value ranging from 0.97P to 1.03P, which results in the minimum repetition time P^obs ₃ ∼ 33 (see the bottom panel in Figure 11). Note that the maximum value of the observed repetition time in the figure is clipped at P^obs ₃ ∼ 150P as ${P}_{3}^{\mathrm{obs}}\to \infty$ for P₃ = 1P when the apparent drift direction changes. The figure clearly shows that we can model the observed nonstable drift characteristic in J1750−3503 with small variations in the actual repetition time when P₃ ∼ 1P.

Zoom In Zoom Out Reset image size

In this section we discuss what is required to change the actual direction of plasma drift with respect to the neutron star rotation.

In the MC model the drift direction depends on the electric potential variation at the polar cap. The spark structure is somehow stable, that is, the distance between sparks does not change significantly while sparks rotate around electric potential extremum. The sparks' location and rotation are forced by the drift of plasma between sparks (see Figure 8 in Szary & van Leeuwen 2017). Note that this idea differs from the idea presented in van Leeuwen & Timokhin (2012) where sparks (plasma columns) are surrounded by vacuum. However, there must be a mechanism that prevents or limits discharges in regions between sparks so they do not result in radio emission in those regions. One of the possible scenarios is a reverse plasma flow, induced by a mismatch between the magnetospheric current distribution and the current injected in the spark-forming regions (see, e.g., Lyubarsky 2012).

There are two main factors that determine the drift direction in the MC model. The first factor is the pulsar geometry and the location of the line of sight with respect to the electric potential extremum within the polar cap. In general, a line of sight passing in between the rotational axis and the electric potential extremum produces a drift direction opposite to that when both the rotational axis and the electric potential extremum are at the same side of the line of sight. Thus, the drift direction can be modified by a change in location of the polar cap or a change in location of the electric potential extremum at the polar cap. A change in polar cap location is possible by a change in structure of the surface magnetic field, while a change in electric potential across the polar cap may be caused by variation in plasma flow from the polar cap region. In the MC model, for an inclined pulsar geometry, a relatively small shift of the order (β/ρ) r_pc may cause drift to change its direction, where r_pc is the polar cap radius. In the case of J1750−3503, assuming purely dipolar configuration of the surface magnetic field, this requires a shift of the order of ∼100 m. Note that if the surface magnetic field is highly nondipolar a very small shift of the order of several meters may cause drift to change its direction. However, for the aligned rotator, such as B0826−34, some other mechanism would have to exist.

The second factor that can influence the drift direction is the type of electric potential extremum at the polar cap. In general, for a pulsar geometry (Ω · B < 0) with a net positive charge at the polar cap, the spark-forming regions, as plasma-starved regions, result in electric potential minimum at the polar cap center. On the other hand, an antipulsar geometry (Ω · B > 0), with a net negative charge at the polar cap, is associated with electric potential maximum. A change in the drift direction requires a significant change in electric potential variation, for example, switching between a minimum and a maximum at the polar cap center, which in principle could be caused by changing the reverse plasma flow in regions between sparks. However, such drastic changes in the electric potential variation on the polar cap would most likely be associated with changes in the sparks structure, which is not visible in a stable pattern of single-pulse emission in J1750−3503 (see Figure 1).

In the LBC model the drift direction is defined by the location of the polar cap with respect to the rotation axis. By introducing a nondipolar surface magnetic field, one can change the size, shape, and location of the polar cap and thus influence the observed drift characteristics (see, e.g., Basu et al. 2020). To change the drift direction, the position of the polar cap with respect to the rotation axis has to change. Thus, it requires a change in configuration of the surface magnetic field. For an aligned rotator, a change in the polar cap position has to be of the order of at least the polar cap size, which for the nondipolar configuration of the surface magnetic field may range from a dozen to several dozen meters (see, e.g., Szary et al. 2017, and references therein). For an inclined rotator, the magnitude of change depends on the initial configuration of the surface magnetic field. If initially an actual polar cap is near the magnetic axis, a change in the polar cap position has to be of the order of at least (α/360°) R. Here α is the inclination angle in degrees and R is the neutron star radius. In the case of J1750−3503 this requires a change in the polar cap position of the order of ∼300–600 m. However, if the initial position of the polar cap is near the rotation axis, then a change in position has to be of the order of the polar cap size, as in an aligned rotator case.

In conclusion, to explain the intrinsic drift direction change, both the MC and LBC models require very rapid significant changes either in electric potential variation across the polar cap or in the position of the polar cap itself.

In Section 3.4 we have shown that the apparent drift direction can be interpreted as an aliasing phenomenon and that we can reproduce the drifting behavior of J1750−3503 with P₃ ∼ P (or less for higher alias order) and small variations of the order of 6%. Both the LBC and MC models have different predictions about the true value of P₃.

In the LBC model, using the simple assumption that LBC is associated with the increasing subpulse phase for successive pulses, the anticorrelation ${P}_{3}={(\dot{E}/\dot{{E}_{0}})}^{-0.6\pm 0.1}$ was found, where E₀ = (2.3 ± 0.2) × 10³² erg s⁻¹ (Basu et al. 2016,2019). That anticorrelation predicts the repetition time of J1750−3503 to be P₃ = (10 ^{+ 6} ₋₃)P. This value seems to be inconsistent with the measured value ${P}_{3}^{\mathrm{obs}}=(43.5\pm 0.4)P$.

It is worth mentioning that if we consider only pulsars with positive drift (according to the definition of drift direction used in this paper) no such anticorrelation can be inferred (see, e.g., the blue circles in the right panel of Figure 1 in Basu et al. 2019). This suggests that some other factors, for instance, the nondipolar structure of the surface magnetic field, may play an important role in determining the actual direction of drift and that further studies are required to verify the suggested correlation. Note that the definition of positive drift in this paper is opposite to the definition used in Basu et al. (2016, 2019) and that there are some errors, for instance, in the caption of Figure 6 in Basu et al. (2016), related to the drift direction definition.

In the original carousel model, on the other hand, the dependence ${P}_{3}\simeq 5.6\,{n}_{\mathrm{sp}}^{-1}\,{B}_{12}\,{P}^{-2}$ was postulated, where n_sp is the number of sparks and B₁₂ is the magnetic field in units 10¹² G (see Equations (33) and (34) in Ruderman & Sutherland 1975). Assuming a dipolar magnetic field at the surface, ${B}_{d}=3.2\times {10}^{19}{(P\dot{P})}^{0.5}\,G$, and the spin-down energy loss $\dot{E}=4{\pi }^{2}{{IP}}^{-3}\dot{P}$, where I ≃ 10⁴⁵ g cm² is the moment of inertia, we can find the dependence P₃ ≃ (5.6/n_sp) ${\left(\dot{E}/{\dot{E}}_{1}\right)}^{0.5}$, where E₁ ≃ 4 × 10³¹ erg s⁻¹. Assuming 13 sparks (see Section 3.4), this results in the repetition time P₃≃ 0.15P, which is consistent with observations assuming alias order n = 7. This high alias order corresponds to the maximum possible drift velocity, that is, the minimum possible value of P₃, and may be unlikely in reality. In the MC model, for instance, the drift velocity depends not on the absolute value, as proposed by Ruderman & Sutherland (1975), but on the variation of the accelerating potential across the polar cap. Thus, in the MC model it is not possible to determine the dependency of P₃ on $\dot{E}$ without knowing the details of sparks formation and existence of regions where they do not arise.

In the case of B0826−34, Gupta et al. (2004) constructed a model that explains its behavior assuming that the actual drift is much faster and we observed only its aliased value. This claim was later called into question by van Leeuwen & Timokhin (2012), who argued that since the observed reversals influence the intrinsic subpulse brightness, the reversals are not just apparent. Furthermore, detailed analysis on a number of other pulsars also concluded their pulse trains were unaliased (see, e.g., van Leeuwen et al. 2003). However, in the MC model, a change in P₃ should be connected with variation in accelerating potential at the polar cap and in principle may have some implications on subpulse brightness. Therefore, we should not exclude the possibility of aliasing in B0826−34 based only on intrinsic subpulse brightness changes.

In the case of J1750−3503 there is no evidence of an intrinsic change of subpulse brightness during reversals (see Section 3.3). However, the average profile characteristics and the separation between subpulses change when drift changes direction (see Section 3.2.2). As pointed out by Gupta et al. (2004), the measured value of P ^m ₂ is a "Doppler-shifted" version of the actual subpulse separation, P ^t ₂, and it depends on the drift rate, P ^t ₂ = P ^m ₂ (1 − D^t/360°), where D ^t is the true value of drift rate. Thus, a varying nature of drift rate is naturally associated with a change in average profile. The drift rate change inferred from the subpulse separation (∼7%) is consistent with a drift rate change required to explain the observed behavior of J1750−3503 by aliasing.

In the recent work of Janagal et al. (2022) it is proposed that the four observed drift modes in J1822−2256 are the result of intrinsic P₃ ∼ P (or less for higher alias order) and a varying number of sparks. It is worth pointing out that although all the drift modes in J1822−2256 have the same direction, in principle different spark configurations (or varying drift rates) in combination with P₃ < 2P could result in opposite drift directions in different modes.