The frequency-dependent behaviour of subpulse drifting: I. Carousel geometry and emission heights of PSR B0031-07

TL;DR

This paper investigates the subpulse drifting behavior of pulsar PSR B0031-07 within the carousel model, explaining its drift modes through changes in spark numbers and aliasing effects, and extends the model to include frequency-dependent emission height estimates.

Contribution

It introduces a unified carousel model explaining all drift modes via spark number changes and aliasing, and develops a new method to estimate emission heights using aberration and retardation effects.

Findings

Drift modes can be explained by a single carousel rotation rate with varying spark numbers.

Confirmed harmonic relation of P3 values across modes.

Estimated emission height differences are less than 2000 km between 185 MHz and 610 MHz.

Abstract

The carousel model of pulsar emission attributes the phenomenon of subpulse drifting to a set of discrete sparks located very near the stellar surface rotating around the magnetic axis. Here, we investigate the subpulse drifting behaviour of PSR B0031-07 in the context of the carousel model. We show that B0031-07's three drift modes (A, B, and C) can be understood in terms of a single carousel rotation rate if the number of sparks is allowed to change by an integral number, and where the different drift rates are due to (first-order) aliasing effects. This also results in harmonically-related values for P 3 (the time it takes a subpulse to reappear at the same pulse phase), which we confirm for B0031-07. A representative solution has [n_A, n_B, n_C] = [15, 14, 13] sparks and a carousel rotation period of P_4 = 16.4 P_1. We also investigate the frequency dependence of B0031-07's subpulse behaviour. We extend the carousel model to include the dual effects of aberration and retardation, including the time it takes the information about the surface spark configuration to travel from the surface up to the emission point. Assuming these effects dominate at B0031-07's emission heights, we derive conservative emission height differences of $\lesssim 2000$ km for mode A and $\lesssim 1000$ km for modes B and C as seen between 185 MHz and 610 MHz. This new method of measuring emission heights is independent of others that involve average profile components or the polarisation position angle curve, and thus provides a potentially strong test of the carousel model.

Tables

Table 1. Table 1: Data sets recorded on MJD 57556

Telescope	Centre frequency	Bandwidth	# pulses	S/N per
	(MHz)	(MHz)		pulse
MWA	185	30.72	4631	14.2
GMRT	610	33.33	8487	9.8

Table 2. Table 2: Values of P 2 subscript 𝑃 2 P_{2} and P 3 subscript 𝑃 3 P_{3} for the drift modes of B0031 − - 07, from McSweeney et al. ( 2017b )

Drift mode	$P_2^{(\circ )}$	${\overline{P}}_3$
A	${18}^{\circ }.9\pm 1.1$	$12.5\pm 0.8$
B	${19}^{\circ }.8\pm 0.5$	$7.0\pm 0.2$
C	${19}^{\circ }.1\pm 2.9$	$4.6\pm 0.3$

Table 3. Table 3: Carousel parameters deduced for B0031 − - 07

$n_A$	$n_B$	$n_C$	$k$	$\sim {\overline{P}}_4$	$\Delta n$
13	12	11	1	14.1	0.89
14	13	12	1	15.2	0.96
15	14	13	1	16.4	1.02
15	16	17	-1	-14.0	-0.87
16	17	18	-1	-14.8	-0.93
17	18	19	-1	-15.7	-0.99

Table 4. Table 4: Carousel rotation rates of selected pulsars

Pulsar	Modes	$n$	$\sim \left|{\overline{P}}_4\right|$		${\overline{P}}_{4\text{,RS}}$	Reference
			$k=0$	$\left|k\right|=1$
B0031$-$07	[A,B,C]	[15,14,13]	[187,98,60]	16.4	4.0	[This work]
B0809$+$74	-	≳15	≳165aa$P_4$ values in the cited texts	≳14	1.6	van Leeuwen et al. (2003)
B0818$-$41	-	20	370aa$P_4$ values in the cited texts	18.3aa$P_4$ values in the cited texts	1.9	Bhattacharyya et al. (2007)
B0826$-$34	-	13	Fluctuates in	13aa$P_4$ values in the cited texts	2.3	Gupta et al. (2004); Esamdin et al. (2005);
			both directions			Bhattacharyya et al. (2008)
B0943$+$10	B	20	43	37aa$P_4$ values in the cited texts	9.3	Deshpande & Rankin (2001)
B1237$+$25	-	18bbDerived for the outer ring of their ‘A’ sequence	∼60	26aa$P_4$ values in the cited texts	3.4	Maan & Deshpande (2014)
B1819$-$22	[A,B,C]	[12,11,10]	[204,80,46]	12.8aa$P_4$ values in the cited texts	2.6	Joshi et al. (2018)
B1857$-$26	-	20	147aa$P_4$ values in the cited texts	18	5.4	Mitra & Rankin (2008)
B1918$+$19	[A,B,C]	[10,9,7]	[61,34,17]	12aa$P_4$ values in the cited texts	7.3	Rankin et al. (2013)

Equations

$$D\equiv \frac{P_2}{{\overline{P}}_3},\overline{D}\equiv \frac{{\overline{P}}_2}{{\overline{P}}_3}\text{and}{D}^{(\circ )}\equiv \frac{P_2^{(\circ )}}{{\overline{P}}_3}.$$

$$\frac{1}{{\overline{P}}_3}=\left|\frac{n}{{\overline{P}}_4}-k\right|,$$

$$\frac{1}{{\overline{P}}_{3B}}=\frac{1}{2}\left(\frac{1}{{\overline{P}}_{3A}}+\frac{1}{{\overline{P}}_{3C}}\right).$$

$$\left|\frac{1}{{\overline{P}}_{3A}}-\frac{1}{{\overline{P}}_{3B}}\right|=\left|\frac{1}{{\overline{P}}_4}\right|.$$

$$\frac{n_A}{\Delta n}=\frac{{\overline{P}}_{3B}(1-k{\overline{P}}_{3A})}{{\overline{P}}_{3B}-{\overline{P}}_{3A}},$$

$$\frac{n_A}{\Delta n}\approx (-1.3\pm 0.2)+(15.9\pm 1.7)k.$$

$${\overline{P}}_{4\text{,RS}}\approx 5.7\times {\left(\frac{P_1}{\text{s}}\right)}^{-3/2}{\left(\frac{\dot{P}}{10^{-15}}\right)}^{1/2}$$

$$\frac{1}{{\overline{P}}_2(\phi )}=\left|\frac{d\theta }{d\phi }\right|.$$

$$\theta (\phi )=n\left(\sigma (\phi )-\frac{\phi }{{\overline{P}}_4}\right)+{\theta }_0,$$

$$\frac{1}{{\overline{P}}_2(\phi )}=n\left|\frac{d\sigma }{d\phi }-\frac{1}{{\overline{P}}_4}\right|.$$

$$\frac{1}{\overline{D}}=\bigg{(}\overline{P}_{3}k\pm 1\bigg{)}\bigg{(}\overline{P}_{4}\frac{d{\sigma}}{d{\varphi}}-1\bigg{)}.$$

$$\frac{{\overline{P}}_{2A}}{{\overline{P}}_{2B}}=\frac{n_B}{n_A},$$

$$\frac{d\sigma }{d\phi }\approx \pm \frac{1}{n{\overline{P}}_2},$$

$${\phi }^{\prime }\approx \phi -2r^{\prime },$$

$${\theta }^{\prime }({\phi }^{\prime })$$

$$=\theta (\phi )+\frac{n}{{\overline{P}}_4}{r}^{\prime }.$$

$${\theta }_{\text{L}}^{\prime }({\phi }_{\text{L}}^{\prime })-{\theta }_{\text{H}}^{\prime }({\phi }_{\text{H}}^{\prime })$$

$$=-\frac{n}{2{\overline{P}}_4}({\phi }_{\text{L}}^{\prime }-{\phi }_{\text{H}}^{\prime }),$$

$$\frac{\Delta {\theta }^{\prime }}{\Delta {\phi }^{\prime }}=-\frac{n}{2{\overline{P}}_4}=-\frac{1}{2}\left(k\pm \frac{1}{{\overline{P}}_3}\right).$$

$$\Delta r^{\prime }\approx -\frac{1}{2}\Delta {\phi }^{\prime },$$

$$\Delta {\phi }_{\text{hor}}^{\prime }\approx \Delta {\phi }^{\prime }+\Delta {\theta }^{\prime }\frac{d{\phi }^{\prime }}{d{\theta }^{\prime }}.$$

$$\Delta {\phi }_{\text{hor}}^{\prime }\approx \Delta {\phi }^{\prime }\left(1-\frac{n{\overline{P}}_2^{\prime }}{2{\overline{P}}_4}\right),$$

$$\Delta r^{\prime }\approx \frac{\Delta {\phi }_{\text{hor}}^{\prime }}{2-k{\overline{P}}_2^{\prime }\pm {\overline{D}}^{\prime }},$$

$$\frac{1}{{\overline{P}}_2^{\prime }}=\frac{d{\theta }^{\prime }}{d{\phi }^{\prime }}\approx \frac{d}{d{\phi }^{\prime }}\left(\theta +\frac{n}{{\overline{P}}_4}{r}^{\prime }\right).$$

$$\frac{d}{d{\phi }^{\prime }}=\frac{d\phi }{d{\phi }^{\prime }}\frac{d}{d\phi }\approx \frac{1}{1-2\frac{d{r}^{\prime }}{d\phi }}\frac{d}{d\phi },$$

$$\frac{1}{{\overline{P}}_2^{\prime }}=\frac{1}{1-2\frac{d{r}^{\prime }}{d\phi }}\left(\frac{1}{{\overline{P}}_2}+\frac{n}{{\overline{P}}_4}\frac{d{r}^{\prime }}{d\phi }\right).$$

$$\frac{1}{{\overline{P}}_2^{\prime }}\approx \frac{1}{{\overline{P}}_2}\left(1+2\frac{d{r}^{\prime }}{d\phi }\right).$$

$$\frac{d{r}^{\prime }}{d{\phi }^{\prime }}=\frac{d\Phi }{d{\phi }^{\prime }}{\nu }^{\eta }$$

$$\frac{1}{{\overline{P}}_2^{\prime }}\approx \frac{1}{{\overline{P}}_2}\left(1+2\frac{d\Phi }{d{\phi }^{\prime }}{\nu }^{\eta }\right).$$

$$\nu_{\text{asymptote}}=\bigg{(}-2\frac{d{\Phi}}{d{\varphi^{\prime}}}\bigg{)}^{-1/\eta},$$

Full text

1 Introduction

Subpulse drifting, first observed by Drake & Craft (1968) not long after pulsars were discovered (Hewish et al., 1968) is a widespread phenomenon affecting more than half of all known pulsars (Weltevrede et al., 2006, 2007). When individual pulses are stacked vertically to form a two dimensional pulsestack, the individual components of pulses, termed subpulses, are observed to undergo a regular modulation pattern in both amplitude and phase, which visually resemble a set of discrete diagonally-oriented burst regions called drift bands. The slope of the drift bands is called the drift rate, and it represents the apparent advance or lag of subpulses per stellar rotation.

Drifting subpulses found an early explanation in the seminal polar gap theory of Ruderman & Sutherland (1975) (hereafter RS), who suggested that the pattern of radio waves that ultimately escape the magnetosphere is the emission signature of a set of discrete, localised pockets of quasi-stable electrical activity called sparks that exist very near the pulsar surface. These sparks are the sites of particles accelerated through the polar gap to energies of $10^{12}\,$V, seeding an avalanche of secondary particles which stream along magnetic field lines at relativistic speeds to produce curvature radiation. Charged particles generated by the sparks partly screen the electric potential along the magnetic field at the surface, causing them to move relative to the polar cap surface about the magnetic axis (the $\mathbf{E}\times\mathbf{B}$ drift), in an arrangement resembling a fairground “carousel”. The emission beam would then also have a spatial structure composed of discrete “beamlets” that reflects the magnetic azimuthal arrangement of sparks on the surface. With each rotation, an Earth-bound observer perceives a different intensity pattern as the line of sight cuts through a slightly rotated carousel, producing the observed drifting behaviour.

The carousel model has been heavily criticised on (at least) two different grounds. Firstly, its reliance on a coherent form of curvature radiation is difficult to justify from known plasma physics (e.g. Melrose, 1992). Secondly, there are a number of secondary phenomena associated with subpulse drifting that present challenges for, or at least require extensions of, the simple geometric picture described above (Edwards, 2006). These include, among others, drift mode switching (Huguenin et al., 1970), bi-drifting (Qiao et al., 2004), nulling (Backer, 1970), and subpulse “memory” across nulls (Lyne & Ashworth, 1983).

Despite these criticisms, the carousel model, along with its implied conal emission beams, remains a popular model for interpretating subpulse drifting, as it is capable of explaining (qualitatively at least) several features of pulsar morphology and behaviour (Rankin, 1983, and subsequent papers in the series). As a result, there have arisen two schools of thought, one which believes the basic model is correct in principle but whose details are complicated and have yet to be worked out on a pulsar-by-pulsar basis, and one which believes the model is fundamentally incorrect. What is lacking is a robust test that can distinguish, once and for all, between these alternatives. A strong candidate for such a test, if it can be found, is one that exploits the model’s unique emission geometry, namely, that although the observed photons are emitted at some height dictated by the magnetic field geometry, the subpulse modulation pattern is dictated by events very close to the surface. This feature of the carousel model makes specific predictions about how the observed subpulse modulation changes as a function of observing frequency.

In this two-part paper series that will focus on B0031$-$07, we investigate the frequency dependence of subpulse drifting. We start by recognising that the basic carousel model predicts that the rotation phases at which subpulses appear is invariant with frequency (Edwards & Stappers, 2002, 2003), a point that is sometimes confused with the expected radius-to-frequency mapping (RFM) behaviour of average profiles (e.g. Smits et al., 2007; Yuen et al., 2016). Nevertheless, frequency-dependent behaviour of subpulses is observed, most notably in the way that the time between successive subpulses (the so-called “secondary period”, $P_{2}$) changes as a function of frequency (Taylor et al., 1975; Bartel et al., 1980). Several ideas have been advanced to explain this within the context of the carousel model, including non-dipolar field geometries (e.g. Davies et al., 1984) and the finite size and shape of beamlets, which present a slightly-shifted cross section to the line of sight as the frequency-dependent beamlets move in magnetic colatitude (Edwards & Stappers, 2003; Bilous, 2018). Edwards & Stappers (2003) also speculate that aberration and retardation (AR) effects must be present, although they argue that the subpulse phase shifts they observe in PSR B0320+39 and PSR B0809+74 imply emission heights that are much larger than the heights expected from other arguments.

In this paper, we develop this idea further to obtain quantitative predictions of how subpulses shift with frequency if AR effects are the dominant cause. We apply the model to the interesting and well-studied case of PSR B0031$-$07 whose $P_{2}$ is also known to have a measurable longitude dependence (McSweeney et al., 2017b) and whose known frequency-dependent subpulse phase shifts are also dependent on which of the pulsar’s three drift modes is present (Huguenin et al., 1970; Vivekanand & Joshi, 1996; McSweeney et al., 2017a). In Sections §3.1 and §3.2 the relationship (initially neglecting AR effects) between the drift bands and the underlying carousel geometry are laid out, with particular emphasis on the degeneracy introduced by the possibility of aliasing effects. We then generalise the model to include AR effects in Section §3.3, calculating their effect on the subpulse phase shift (§3.3.1) and on $P_{2}$ (§3.3.2). A discussion of the results and our final conclusions are presented in Sections §4 and §5.

This is the first of a two-part series, both of which deal with the relationship between frequency-dependent effects and subpulse drifting in B0031$-$07. This paper deals exclusively with effects that can be observed in Stokes I observations. In the subsequent paper, we investigate how the subpulse phase (which encodes information about the magnetic azimuth) is related to the polarisation position angle and use this relationship to place further constraints on the carousel and viewing geometries.

1.1 Definitions and Notation

The rotation axis ($\vec{\omega}$), the magnetic axis ($\vec{\mu}$), and the LoS define three points of a spherical triangle which relate the rotation phase ($\varphi$), the viewing geometry ($\alpha$, $\zeta$), the half-beam opening angle ($\Gamma$), the PPA ($\psi$), and the magnetic azimuth ($\sigma$). In this paper, $\alpha$ and $\zeta$ are always relative to the rotation axis defined by the right hand rule, so that $0^{\circ}\leq\alpha<180^{\circ}$ and $0^{\circ}\leq\zeta<180^{\circ}$. These angles are illustrated in Fig. 1. An asterisk is used to denote the supplementary angles $\alpha^{\ast}=180^{\circ}-\alpha$ and $\zeta^{\ast}=180^{\circ}-\zeta$, i.e. measured from the opposite rotation axis. The impact angle is defined as $\beta=\zeta-\alpha$, with $\beta^{\ast}=\zeta^{\ast}-\alpha^{\ast}=-\beta$.

As is common in the literature, we use $P_{1}$ for rotation period of the pulsar, $P_{2}$ for the subpulse separation, $P_{3}$ for the time it takes a subpulse to reappear at the same rotation longitude, and $P_{4}$ for the rotation period of the carousel. Conventionally, $P_{2}$ is often expressed in degrees of rotation, and $P_{3}$ in units of $P_{1}$, but in order to keep the notation consistent throughout this work, we reserve the unadorned $P_{\text{x}}$ to have units of time. To denote quantities normalised by $P_{1}$, we introduce the overline notation, e.g., $\overline{P}_{3}\equiv P_{3}/P_{1}$. In addition, we define $P_{2}^{(\circ)}\equiv\overline{P}_{2}\times 360^{\circ}$ to express periods in degrees of rotation.

The drift rate, $D$, is traditionally expressed in units of degrees per pulse. To keep this consistent with the foregoing notation, we define

[TABLE]

2 Observations

We have taken two simultaneous observations of B0031$-$07, one using the Murchison Widefield Array (MWA; Tingay et al., 2013) centred at $185\,$MHz, and the other using the Giant Metrewave Radio Telescope (GMRT; Swarup et al., 1991), centred at $610\,$MHz. The details are summarised in Table 1.

The MWA observation was processed in the same way as McSweeney et al. (2017b). The MWA voltages were recorded using the Voltage Capture System (VCS; Tremblay et al., 2015), and a tied-array beam was formed (Ord et al., 2019) in the direction of the pulsar. The calibration solution was obtained using a separate, dedicated observation of the calibrator source 3C444. During the calibration process, it was deemed necessary to exclude five MWA tiles (out of 128) from the coherent sum, due to corrupted data. The resulting data set had a time resolution of $100\,\mu$s and a frequency resolution of $10\,$kHz. At the time of the observation (June 2016), only total intensity data products were available, but the recent verification of MWA polarimetry (Xue et al., 2019) has made it possible to reprocess the original recorded voltages to produce a data product with full Stokes parameters.

The GMRT observations were made with the 13 central antennas in the phased array mode. However, only the total intensity was recorded. The output time resolution was $123\,\mu$s and the frequency resolution was $65\,$kHz.

Both MWA and GMRT data sets were processed with DSPSR (van Straten & Bailes, 2011) and PSRCHIVE (Hotan et al., 2004; van Straten et al., 2011) to produce single-pulse archives. A small percentage of frequency channels contaminated by local radio frequency interference (RFI) were removed, and impulsive RFI was identified by eye on a pulse-by-pulse basis; pulses that were contaminated by RFI were replaced by zeros. The resulting Stokes I average profiles are shown in Fig. 2.

The analysis in this paper will naturally benefit from knowledge of the exact clock offset between the two telescopes, which unfortunately was unavailable at the time of recording due to the fact that the MWA had not yet embarked on any timing programs that could verify the VCS time stamps. The correspondence between pulses was therefore found by aligning the null sequences, which are known in the case of B0031$-$07 to be a broadband phenomenon (Smits et al., 2007). The total number of pulses observed simultaneously at both telescopes was $4780$. A section of the pulsestacks, showcasing all three drift modes, is shown in Fig. 3.

3 Analysis & Results

3.1 Using $P_{3}$ to constrain the carousel geometry

In the simplest version of the carousel model, in which the surface sparks are equally spaced in magnetic azimuth, the behaviour of the drift bands is determined by the number of sparks in the carousel, $n$, and the carousel rotation rate, $P_{4}$. If the carousel rotation is constant, the observed $P_{3}$ at any given longitude will also be constant. This is true regardless of either the radial distance of the spark from the magnetic axis, or the sparks’ radial extent. Without aliasing, $P_{3}$ is identical to $P_{4}/n$, but with aliasing,

[TABLE]

where $k=[n/\overline{P}_{4}]$ is the aliasing order111Under this definition of the aliasing order, the sign of $k$ is always the same as the sign of $P_{4}$, which in this work is chosen to be positive when the carousel is rotating anticlockwise when viewed from “above” the magnetic axis. The phrase “first order aliasing” can be taken to mean $k=\pm 1$., and where the square brackets denote rounding to the nearest integer. Eq. (2) is the same as that derived by Rankin et al. (2013), but generalised for arbitrary aliasing orders.

The fact that $P_{3}$ is different for each of B0031$-$07’s three drift modes can therefore be explained by a change in $n$, $P_{4}$, or both. In general, it is difficult to find the correct aliasing order $k$, because solutions for Eq. (2) can always be found for any $k$ due to the free choice of $P_{4}$ for each mode. However, the carousel rotational speed is thought to be set by the magnetic and electric fields near the surface, which cannot easily change their magnitudes and/or directions on the time scales on which drift mode changes are observed to take place ($<1$ rotation). Therefore, it is an interesting question to ask whether the behaviour of the different modes might be brought about by only a change in $n$, and keeping $P_{4}$ constant, an idea first introduced by Rankin et al. (2013) for PSR B1918+19. In this case, it is immediately clear that $k$ cannot be [math] for all three drift modes, since otherwise the drift rate would be the same at all times. If we assume that $k$ is the same for all three drift modes, we find that the quantities $1/\overline{P}_{3A}$, $1/\overline{P}_{3B}$, and $1/\overline{P}_{3C}$ (for modes A, B, and C, respectively) are members of an arithmetic sequence whose $n$th term corresponds to a carousel with $n$ sparks. Furthermore, if $n_{A}$, $n_{B}$, and $n_{C}$ are themselves an arithmetic sequence, then the respective values of $\overline{P}_{3}$ will be harmonically related, as originally suggested by Wright & Fowler (1981). It also immediately follows that

[TABLE]

which can be readily verified in the case of B0031$-$07. In the special case that the number of sparks differs only by one, then the terms of the arithmetic sequence are consecutive, and the difference between any adjacent pair yields the carousel rotation rate directly, e.g.

[TABLE]

Then $n$ can be found by substituting the derived value of $\overline{P}_{4}$ back into Eq. (2).

An alternative approach rearranges Eq. (2) to leverage the expectation that the number of sparks changes by an integer value. Continuing under the assumption that the carousel rotation does not change between drift modes, $\overline{P}_{4}$ can be eliminated to yield

[TABLE]

where $\Delta n=n_{A}-n_{B}$. Using measured values for $\overline{P}_{3A}$ and $\overline{P}_{3B}$ (from McSweeney et al., 2017b, summarised in Table 2), this becomes

[TABLE]

Plugging in candidate values for $n_{A}$ and $k$, we find near-integral solutions for $\Delta n$, with the first order solutions given in Table 3. As expected, no solutions were found for $k=0$. The same procedure was repeated using modes B and C in Eq. (5) instead of A and B, with the same results.

Eq. (6) shows that if a change in $n$ is indeed the cause of different drift modes, then the predicted number of sparks must increase rapidly for higher aliasing orders, and more so for higher values of $\Delta n$. If $\Delta n$ is kept constant, then $n$ increases approximately linearly with $k$, and therefore $\overline{P}_{4}$ must stay relatively constant (since $k=[n/\overline{P}_{4}]$), but will be much larger (i.e. the carousel rotates much slower) for higher values of $\Delta n$. Thus, the approximate solution for $\overline{P}_{4}$ represents the fastest possible carousel rotation speeds consistent with the observed $\overline{P}_{3}$ values, without significant difference for higher aliasing orders. By contrast, the theoretical value of $\overline{P}_{4}$ predicted by the RS model is

[TABLE]

which is $\overline{P}_{4\text{,RS}}\approx 4$ for B0031$-$07. Thus, we argue that a modest carousel rotation period ($14\lesssim|\overline{P}_{4}|\lesssim 17$) and a modest number of sparks ($13\leq n_{A}\leq 17$) is the most plausible, since it is much closer to the above theoretical value than a no-aliasing solution (for instance, $\overline{P}_{4}=n\overline{P}_{3}=9\times 12.5=112.5$ for mode A in Smits et al. 2005). These carousel properties are similar to those found by Rankin et al. (2013) for PSR B1918+19.

On the possible harmonicity of the drift modes, we note that the finding of Vivekanand & Joshi (1996) that they are not harmonically related appears to be based on the argument that the drift rates are not harmonically related. However, this is not the same as the claim that the values of $P_{3}$ are harmonically related. Indeed, since $P_{2A}\approx P_{2B}\approx P_{2C}$, the harmonicity of $P_{3}$ implies that the drift rates should be linearly related, which is in fact consistent with their measurements.

3.2 Using $P_{2}$ to constrain the viewing geometry

The preceding section dealt solely with $P_{3}$, a quantity which is independent of pulse longitude, as evident in Eq. (2). $P_{2}$, on the other hand, encodes information about the angular distance between successive sparks, and is generally a longitude-dependent quantity. We now proceed to investigate how $P_{2}$ may be used to constrain the carousel parameters and the viewing geometry.

$P_{2}$ can be measured by correlation methods (e.g. Smits et al., 2005) or methods using fluctuation spectra (Backer, 1971; Wolszczan et al., 1981). For our purposes, the latter are the most suitable for obtaining measurements of $P_{2}$ as a function of longitude. This involves computing the longitude-resolved fluctuation spectrum (LRFS), which is done by applying the Discrete Fourier Transform to each longitude bin within a given drift sequence. If the drift band separation remains constant throughout the sequence, a narrow, strong feature will appear at all longitudes in the frequency bin corresponding to the value of $P_{3}$. The relative phases of the feature, which we denote by $\theta(\varphi)$, constitute the so-called “subpulse phase track”, which encodes information about $P_{2}$:

[TABLE]

The subpulse phase that is visible at a given rotation phase involves the computation of two quantities: (1) the magnetic azimuth, $\sigma$, of the field lines coplanar with the LoS, and (2) the spark activity at the footpoint (i.e. where the magnetic field line intersects the stellar surface). The consideration of (1) is what Yuen et al. (2016) mean by the “motion of the visible point.” For a carousel of $n$ sparks rotating at rate $1/\overline{P}_{4}$, the subpulse phase is the difference between the sampled magnetic azimuth and the amount that the carousel has rotated in the same time:

[TABLE]

where the dependence on pulse phase has been made explicit, and $\theta_{0}$ is the subpulse phase when $\varphi=0$. An example of a mode B sequence, observed simultaneously with the MWA and the GMRT, with their respective profiles, LRFS, pulse energies, and subpulse phase tracks, is shown in Fig. 4.

To recover the expression for $\overline{P}_{2}$, we differentiate with respect to the phase, yielding222One must be careful when comparing this equation with similar-looking equations in Yuen et al. (2016), as both the definitions and notations of the terms involved, as well as the derivation itself, differ.

[TABLE]

By multiplying through by $\overline{P}_{3}$ and by using Eq. (2) to eliminate $n$, the drift rate can be equivalently expressed as

[TABLE]

In this form, it is easily seen that if $k=0$, the dependence on $P_{3}$ vanishes and a constant carousel rotation speed implies the same drift rate for all three drift modes. Furthermore, any observed curvature of the driftbands is necessarily inherited from the $d\sigma/d\varphi$ term unless the angular speed of the sparks is a function of magnetic azimuth.

If the drift mode changes are, in fact, characterised by only a change in the number of sparks, then Eq. (10) predicts only small changes in $\overline{P}_{2}$ for sufficiently large values of $n$. This is because $\overline{P}_{4}$ is assumed constant, and both $\sigma$ and $\frac{d{\sigma}}{d{\varphi}}$ are pure functions of the viewing geometry. Therefore, the fractional change in $\overline{P}_{2}$ should be equal to the fractional change in $n$, or

[TABLE]

and similarly for modes A and C, etc. The smallest integral values of $n$ that are consistent with the reported values of $P_{2}$ in McSweeney et al. (2017b) are approximately $n\gtrsim 6$. This is a weaker constraint than the constraint derived using $P_{3}$ because of the relatively large measurement errors on $P_{2}$ (obtained using correlation methods).

The relative magnitudes of the two terms in the brackets of Eq. (10) determine the overall drift band behaviour. For a sufficiently slowly rotating carousel, the value of $P_{2}$ is dominated by $d\sigma/d\varphi$. If the carousel rotates in the opposite sense to the motion of the visible point in the pulse window (i.e. $P_{4}$ is negative), then increasing the carousel speed causes $P_{2}$ to decrease. However, if the carousel rotates in the same direction as the motion of the visible point, $P_{2}$ will initially increase until the two terms balance and the subpulse drift appears stationary. As the carousel rotates faster still, the sparks overtake the visible point and the subpulse drift resumes, but now with each successive subpulse coming from the trailing spark.

Regardless of the carousel configuration, $d\sigma/d\varphi$ is in general a function of rotation phase, and thus so is $P_{2}$. It follows that the drift rate can never be stationary at all phases simultaneously. In the presumably rare case that the value of $1/\overline{P}_{4}$ falls between the maximum and minimum value of $d\sigma/d\varphi$, one would see the subpulse drift adopt opposite directions at distinct rotation phases, a phenomenon that has been observed in a few pulsars and known as “bi-drifting” (Wright & Weltevrede, 2017, and references therein). Thus, Eq. (10) potentially provides a natural framework for understanding bi-drifting.

In the case of B0031$-$07, we can estimate the relative contributions of the two competing terms in Eq. (10) by considering what parameter values give rise to the measured $P_{2}$ value. Taking, for example, the B mode of the second line in Table 3 ($n=13$, $\overline{P}_{4}=15.2$), a value of $d\sigma/d\varphi=1.46$ is required to reproduce $\overline{P}_{2}^{\circ}=19.8^{\circ}$. Thus, even though the carousel is rotating fast enough to produce aliasing $k=1$, the motion of the visible point is far more dominant, since $d\sigma/d\varphi\gg 1/\overline{P}_{4}$. In this case, we can then write

[TABLE]

which becomes an exact equality in the stationary carousel limit. For the same viewing geometry (i.e. keeping $d\sigma/d\varphi$ fixed), the contribution by the carousel rotation only becomes appreciable when either the aliasing order increases to $k\gtrsim 5$ or the number of sparks drastically increases (or some combination of both). A similar conclusion is reached for all the entries in Table 3.

If we assume a small $k$ for B0031$-$07, we may draw another conclusion from the observation that in the pulses immediately following the onset of a drift sequence after a null, the drift rate appears to undergo a decay until it reaches a nominally stable value (McSweeney et al., 2017b). During this decay period, both $P_{2}$ and the subpulse width appear visually exaggerated, suggesting that the relative speed of the visible point and the sparks is less. However, when the motion of the visible point dominates, there are only two possibilities. Either (1) the directional sense of both the carousel drift and the motion of the visible point is the same, in which case the relaxation of the drift rate implies a slowing down of the carousel rotation, or (2) the direction senses are different, in which case the carousel must be speeding up. Given that the drift rate of other pulsars which are thought to be free of aliasing effects is observed to increase after a null (e.g. Lyne & Ashworth, 1983), we therefore favour the second possibility. Note that this does not necessarily imply that the carousel and the pulsar are rotating in the same directional sense, since the sign of $d\sigma/d\varphi$ depends, in part, on the sign of $\beta$, which is still unknown.

3.3 Generalisation for AR effects

The $P_{2}$ model derived in the previous section has no frequency dependence whatsoever. However, $P_{2}$ has long been known to be a function of frequency for pulsars generally (Taylor et al., 1975) and for B0031$-$07 in particular (Bartel et al., 1980). More recently, McSweeney et al. (2017a) showed that the frequency-dependent shift in subpulse phase also depends on which drift mode is present.

One might expect that the dependence of $P_{2}$ on observing frequency is a simple side-effect of RFM, where the “space” between the drift bands grows in proportion to the opening angle of the emission cone (e.g. Smits et al., 2007; Yuen et al., 2016). However, it cannot be so, as discussed at length in Edwards & Stappers (2002, 2003). While it is true that the line of sight would cut through the emission cone at different phases depending on the frequency, the image of a single spark projected onto the pulsar’s sky would not change relative to the line of sight traverse. Observationally, this implies that different parts of drift bands are illuminated when viewed at different frequencies, not that the drift bands themselves shift in phase.

Despite the preceding argument, a radius-to-frequency mapping remains a necessary (if not a sufficient) condition to explain the $P_{2}$ dependence on frequency. For example, Bilous (2018) explains the drift band shifts observed in PSR B0943+10 in terms of the offset phase of the centroid of the line of sight traverse through a beamlet of finite size and shape, as the beamlet changes its magnetic colatitude with frequency.

In this section, we consider the possibility that the frequency-dependence of the subpulse phase shift and of $P_{2}$ is due to aberration and retardation (AR) effects, as suggested by Edwards & Stappers (2003). If AR are the dominant effects, such a generalised model can be used to extract information about the emission heights from multi-frequency measurements of the subpulse phase.

3.3.1 Derivation of the subpulse phase shift

Dyks et al. (2004) and others showed that pulsar emission will be shifted earlier in phase according to

[TABLE]

where $\varphi^{\prime}$ is the observed phase, $r^{\prime}=r/r_{L}$, and $r_{L}=cP_{1}/(2\pi)$ is the light cylinder radius. One $r^{\prime}$ comes from the aberration of the emission angle from the tangent to the field line, and the other comes from retardation, or the finite time it takes the light to travel from the emission point to the observer.

Although Eq. (14) applies to all emission, a correct determination of the observed subpulse modulation must also account for the finite time it takes the information about the subpulse phase to travel from the surface to the emission height. To first order, the trajectory of relativistic particles climbing from the surface to the emission point along a corotating field line is simply $r$, with a corresponding travel time equal to a phase rotation of $r^{\prime}$. This means that the observed subpulse phase does not represent the spark activity at the footpoint of the field line at the time of emission, say, at phase $\varphi$, but rather the spark activity at the time when the phase was $\varphi-r^{\prime}$. Thus, emission that comes from a given magnetic azimuth will be observed at phase $\varphi^{\prime}\approx\varphi-2r^{\prime}$ and have subpulse phase

[TABLE]

Let L and H represent two field lines that have the same magnetic azimuth, $\sigma$, but whose footpoints have different magnetic colatitudes. In the absence of AR effects, the two lines come into view at the same phase, $\varphi$, at which moment the geometry is represented by the same fundamental triangle. Let L represent the line whose footpoint has the smaller colatitude. For a dipolar field, one expects from geometric arguments (neglecting rotation effects on curvature, as discussed in Thomas & Gangadhara 2007 and Thomas et al. 2010) that the lower frequency emission will be observed from L and the higher frequency emission from H, and that the emission occurs at a greater altitude for L than for H. From Eq. (15),

[TABLE]

where we have used Eq. (14) in the second equation. The corresponding point on the subpulse phase track will appear to have moved vertically by $\Delta\theta^{\prime}=\theta^{\prime}_{\text{L}}(\varphi^{\prime}_{\text{L}})-\theta^{\prime}_{\text{H}}(\varphi^{\prime}_{\text{H}})$ and horizontally by $\Delta\varphi^{\prime}=\varphi^{\prime}_{\text{L}}-\varphi^{\prime}_{\text{H}}$ when viewed simultaneously at two distinct frequencies. Curiously, the slope of this apparent motion is constant:

[TABLE]

The shift of a whole subpulse phase track is schematically illustrated in Fig. 5.

If $k$ is known, the two phase tracks can be directly compared, and the phase difference gives a direct measure of the difference of the emission heights at the two frequencies:

[TABLE]

in accordance with Eq. (14). Note that although the direction of the subpulse phase track shift is fixed by the carousel geometry, the magnitude of the shift depends on $r^{\prime}$, which may well be a strong function of pulse phase (e.g. Gangadhara, 2004; Thomas et al., 2010). This can introduce further distortions in the shape of the subpulse phase track as it appears at different frequencies.

The phase shift of individual subpulses can be estimated by Taylor expanding the rotation phase shift as a function of the corresponding subpulse phase shift and the slope of the subpulse phase track:

[TABLE]

By applying Eq. (16) and identifying $\overline{P}^{\prime}_{2}=d\varphi^{\prime}/d\theta^{\prime}$ as the AR-corrected subpulse separation actually measured by the observer, we find

[TABLE]

which by virtue of Eq. (18) can be rearranged to yield

[TABLE]

which directly relates the phase shift of individual subpulses to the difference of emission heights at two simultaneously observed frequencies.

3.3.2 Derivation of $P_{2}$

The preceding analysis can be recast explicitly in terms of $P_{2}$ rather than the subpulse drift phase. It follows from Eq. (15) that the measured value of $\overline{P}_{2}^{\prime}$ (where, as before, the primes on $\overline{P}_{2}^{\prime}$ and $\varphi^{\prime}$ denote values measured in the observer’s AR-corrected frame) is

[TABLE]

In order to relate this to the $\overline{P}_{2}$ of Eq. (10), we may apply the chain rule, noting that

[TABLE]

where we have taken the derivative of Eq. (14) with respect to $\varphi$. Although it is possible to express these equations in terms of $dr^{\prime}/d\varphi^{\prime}$, which is the rate of change of emission height with AR-corrected longitude as it would appear in the observer’s frame, $dr^{\prime}/d\varphi$ is easier to relate to the emission geometry. Then Eq. (22) becomes

[TABLE]

Two simplifying assumptions can be made. First, as before, we expect the $n/\overline{P}_{4}$ term to be small compared to $1/\overline{P}_{2}$. Second, $dr^{\prime}/d\varphi$ itself must be small on average, since the maximum height is not expected to be more than a few percent of the light cylinder radius. If $\beta$ is very small, then $r^{\prime}$ may change rapidly near the fiducial point, but otherwise we assume $dr^{\prime}/d\varphi\ll 1$. Thus, neglecting the $n/\overline{P}_{4}$ term and Taylor expanding the fraction,

[TABLE]

Thus, the observed subpulse separation includes both a fixed, frequency-independent term equal to the subpulse separation that would be observered in the absence of AR effects, and a frequency-dependent term that depends on the rate of change of the observed emission height with rotation phase.

In B0031$-$07, $\overline{P}_{2}^{\prime}$ follows the same trend as the separation of average components, i.e. it becomes smaller at higher frequencies. Thus, the frequency-dependent term must be positive. Since the smallest emission height is expected to be observed at the fiducial point, this supports the idea, suggested by McSweeney et al. (2017b), that the fiducial point lies on the leading side of the pulse window.

Eq. (25) is almost in a form that can be used to fit historical measurements of $P_{2}$ against observing frequency. If the RFM is independent of magnetic azimuth, the emission height can be written in the form $r^{\prime}=\Phi(\varphi)\nu^{\eta}$, where $\Phi(\varphi)$ encapsulates the fact that the magnetic colatitudes of the footpoints of the observed field lines is a function of phase, as illustrated in the geometry of Gangadhara & Gupta (2001) and Gangadhara (2004). It then follows that, for a given fixed observing frequency $\nu$,

[TABLE]

and

[TABLE]

3.4 Application of the model to B0031$-$07

Frequency-dependent subpulse behaviour can be thought of in terms of an absolute phase shift (a bulk shift towards earlier or later phases) and a changing $P_{2}$ (a compression or expansion of the subpulse modulation pattern). However, measuring the absolute phase shift for our dual-frequency observation of B0031$-$07 proved problematic because the exact time offset between the clocks at the two telescopes was unknown. The average profiles for each drift mode are known to evolve significantly with frequency, so they cannot be used to align the phases reliably333A wide-band observation using a single telescope would obviate this problem, unless the profile evolution interferes with the correct determination of DM.. Alternatives, such as methods involving the alignment of single pulses (e.g. via cross-correlating the time series), would mask any average subpulse phase shift that does exist between the two frequencies. Therefore, the best that can be done with our present data sets is to show how the estimated emission height difference depends on an assumed clock offset between the two telescopes. Using $k=1$ and the measured average values of $\overline{P}_{2}$ and $\overline{P}_{3}$ for the three drift modes, the relative emission height differences between the two frequencies are plotted for each mode in Fig. 6. The measured value of $\varphi^{\prime}_{\text{hor}}$ was derived by cross correlating each pulse at the two frequencies, and average all the resulting correlation curves for pulses belonging to the same drift mode. The resulting average correlation curves were smoothed with a Savitzky-Golay filter (Savitzky & Golay, 1964) using a fifth order polynomial and then using cubic-spline interpolation to identify the location of the primary peak.

Although the MWA-GMRT clock offset is unknown, we mark two offsets of interest in Fig. 6. The first corresponds to assumption that the peaks of the profiles at the two respective profiles correspond to the same rotation phase. This offset defines the zero point of the $x$-axis used in the figure. The second offset uses the peak of a particularly bright, narrow burst that appears at both frequencies (pulse number 579 in the MWA data set, equivalent to pulse number 149 in the GMRT data set) to align the clocks. The possibility that the two pulses occur at the same “true” rotation phase is founded on the identification of the bright emission as a candidate “giant pulse”, which are known to come from very narrow phase ranges at some frequencies (Kuzmin & Ershov, 2004; Tuoheti et al., 2011; Nizamdin et al., 2011). If a broadband burst comes from a sufficiently small volume of the magnetosphere, as would be the case for a highly localised surface spark event that generates a large range of particle energies, then the phase shift due to AR effects would apply to all frequencies equally, and the burst would appear at a fixed longitude at all frequencies. However, since Fig. 6 shows that this assumption leads to negative height differences in this model (the vertical dashed line at $-4.2^{\circ}$), we conclude that the bright emission must also come from a range of heights, thus making the derived clock offset erroneous. We note the possibility that the GMRT bright pulse is actually a burst of RFI, as it resembles other instances of off-pulse RFI found in the data set, but we do not investigate this further.

Fig. 6 shows that the height differences for all three modes are positive when the MWA-GMRT phase offset $\gtrsim-0.^{\circ}5$. In this region, the mode with the greatest $\Delta r$ is mode A, with the height differences in modes B and C being smaller but similar to each other. For phase offsets less than $5^{\circ}$, $\Delta r\lesssim 2000\,$km for mode A and $1000\,$km for modes B and C (all within a few percent of the light cylinder radius). In contrast, the height differences found by Smits et al. (2007) were $\sim 10\,$km across a much wider frequency range ($157\,$MHz to $4850\,$MHz); however, their model was based (in part) on the idea that $P_{2}$ is a direct manifestation of RFM, which, as we have argued above, is an incorrect interpretation.

The unknown clock offset does not affect measurements of $P_{2}$, which is related only to the relative phase shift of subpulses at different longitudes. Fig. 7 shows a fit for historical measurements of $P_{2}$ during mode B drift sequences (mode B is the most common mode with the most reliable data, and has been observed at a much wider range of frequencies than either mode A or mode C), over a wide range of frequencies. There is a clear general trend that $P_{2}$ gets smaller with increasing frequency, as previous authors have noted. As indicated in the figure, some of the early, low-frequency measurements of $P_{2}$ relied on single polarisation telescope feeds. It is unclear to what extent this would affect the measurements, given the known complexity of B0031$-$07’s subpulse polarisation (Taylor et al., 1975; Gould & Lyne, 1998) and similar examples of $P_{2}$ ambiguity in other pulsars (e.g. Rankin et al., 2005). However, even if these particular observations were omitted, we find that there is still strong support for a frequency-dependent (i.e. non-constant) $P_{2}$, especially if one accepts the data from Smits et al. (2007), whose method for measuring $P_{2}$ was identically used for (sometimes simultaneous) observations across a broad range of frequencies, albeit using three different telescopes. We thus proceed on the cautious assumption that all the historical measurements presented in Fig. 7 are valid. Using Eq. (27), the fitted parameters are $\overline{P}_{2}^{(\circ)}=23.^{\circ}6\pm 3.^{\circ}8$, $d\Phi/d\varphi^{\prime}=0.3\pm 0.1$, and $\eta=0.70\pm 0.34$, which is also plotted on the figure for comparison. A discussion of this result and its implications is deferred until Section §4.2.

4 Discussion

4.1 Carousel geometry

We have attempted to find plausible, quantitative explanations for two interesting features of B0031$-$07’s subpulse modulation pattern: (1) the changing $P_{3}$ and drift rate associated with its three distinct drift modes (A, B, and C), and (2) the dependence of subpulse arrival phases on observing frequency and on drift mode. The interpretation of both of these effects depend on the carousel model, assumed in this analysis, which connects the observed modulation pattern of subpulses to a set of sparks near the pulsar’s surface rotating around the magnetic axis. Measuring the number of sparks, $n$, and the carousel rotation period, $P_{4}$, is an important step in verifying the predictions of the carousel model. However, these measurements are difficult to make, as there are often degenerate sets of parameters that produce similar looking drift bands in the pulse window. We have argued that the most likely carousel configuration is one in which the subpulses are aliased ($k=\pm 1$), with the number of sparks being $n\approx 15$, and where $n_{A}$, $n_{B}$, and $n_{C}$ differ by $1$ in turn.

This solution explains a number of otherwise puzzling features of B0031$-$07’s drift modes. Firstly, the relationship between the modes’ respective $P_{3}$ values has long been suspected to be a harmonic one. We have shown that this is a mathematical consequence of the number of sparks differing by a constant amount (in this case, one). Also, the similarity of $P_{2}$ in all three modes favours solutions with a relatively large number of sparks, as per Eq. (12), which is incompatible with a no-aliasing solution by virtue of Eq. (6). This solution has the added benefit that the different drift rates of the three modes can be understood in terms of a changing number of sparks rather than a change in the carousel drift rate itself (Rankin et al., 2013). This is more consistent with the original claim of RS that the carousel speed is a (constant) function of the electric and magnetic field strengths at the surface, which are not expected to change rapidly on small time scales.

The absolute magnitude of the carousel rotation period is also an important prediction of the RS and competing models. The solutions presented in Table 3 all predict a similar carousel rotation period; taking one example ($[n_{A},n_{B},n_{C}]=[15,14,13]$), the carousel rotation period is $\overline{P}_{4}=\pm 16.4$. The time (per rotation period) it takes one spark to reach the location of its predecessor is of order $\overline{P}_{4}/n\approx 1$, as expected for a solution with first-order aliasing. This is much smaller than the time that would be inferred directly from $\overline{P}_{3}$ if aliasing were assumed to be absent ($12.5$, $7.0$, $4.6$ for modes A, B, and C respectively).

In Table 4, we compare our measurement of $P_{4}$ (for a particular carousel configuration from Table 3) alongside measurements for other pulsars. The carousel rotation period predicted from RS, $\overline{P}_{4\text{,RS}}$, is calculated from Eq. (7). We have included both the unaliased and (first order) aliased values of $P_{4}$ for each pulsar in the table except B0826$-$34, whose drift bands frequently change directions, and whose unaliased $P_{4}$ value is therefore ambiguous. The preferred $P_{4}$ (according to the cited work) is indicated in the table notes, with the other value being extrapolated assuming the same number of sparks. We have not indicated whether $k=1$ or $k=-1$, but this can be worked out using the table values and Eq. (2).

In almost all cases, the unaliased carousel rotation period is much larger than the value predicted by RS, a point that is noted by all authors who advance the unaliased model as the preferred one. It is commonly recognised that aliasing could be a way to account for (or at least mitigate) this discrepancy, except in the case of B0809$+$74, in which the “memory” of subpulse phases across null sequences was used to show that only the unaliased model is consistent with the data. Even if aliasing were possible for this and other pulsars, there still remains a large discrepancy between the measured $\overline{P}_{4}$ and the RS value, most significantly for B0809$+$74, B0818$-$14, and B1237$+$25, where there is approximately an order of magnitude difference. It is possible that in some of these cases, an even higher aliasing order may be present.

Alternative theories of spark circulation obviate the need for invoking higher order aliasing, such as the spark model of Gil & Sendyk (2000), which generally predicts much longer carousel rotation periods than RS. They analyse B0031$-$07, concluding that its carousel consists of five sparks in a single annulus with a radius of $s=0.7$ of the polar cap distance, rotating at $\overline{P}_{4}\approx 34$. While this reproduces an unaliased value of $P_{3}$ that agrees well with the observed value for $P_{3}$, they do not discuss how their model may be applied to modes A and C. However, we note that their theory predicts a maximum allowed number of sparks of $n_{\text{max}}=16$, which is consistent with the $k=+1$ solutions in Table 3.

A final point to note is that, in many cases, the first order aliasing solution predicts a value of $\overline{P}_{4}$ that is comparable to the number of sparks, and hence $n/|\overline{P}_{4}|\sim 1$. This is not surprising, as first order aliasing implies $0.5\leq|k|=[n/|\overline{P}_{4}|]\leq 1.5$.

4.2 Frequency-dependent effects

Armed with a plausible solution for the configuration and rotation period of the carousel, we have attempted to make sense of the frequency-dependence of its subpulse modulation pattern. It is important to remember that the carousel model (ignoring AR and finite-spark-size effects) predicts that any given subpulse should appear at the same rotation phase at any observing frequency, regardless of the location (or movement) of their parent sparks in the polar cap, and regardless of the emission height. Clearly, however, the subpulses of B0031$-$07 do exhibit frequency- and mode-dependent phase shifts, as has been demonstrated many times in the literature, and confirmed by our own observations.

The qualitative behaviour of the subpulses can be summarised as follows. In any given drift mode, the horizontal separation of the drift bands (i.e. $P_{2}$) decreases with increasing frequency (however, mode A may behave more erratically at higher frequencies, see Smits et al., 2007) by roughly the same amount. There is also an absolute, frequency-dependent phase shift of the whole subpulse pattern which is different for the three modes.

4.2.1 Absolute phase shift and implications for emission heights

We have derived a model for the frequency-dependence of subpulse phases due to AR effects. Regardless of whether this turns out to be the correct explanation for B0031$-$07, the derived formulas quantify the subpulse phase shift that must be present due to AR effects. Importantly, we have included in this model the time it takes information about the spark modulation to travel from the surface to the emission height, a point which (to the best of our knowledge) has not been considered in previous studies of the frequency-dependence of subpulses. This makes the model a potentially powerful tool for testing the carousel model in cases where AR effects are indeed dominant, and where the emission heights can be measured by other (independent) means.

Our attempt to apply this model to B0031$-$07, however, encountered several difficulties. Firstly, the exact time offset between the data obtained at the two different telescope facilities is unknown, so that the absolute phase shifts of the subpulse modulation pattern in each drift mode could not be confidently measured. This analysis is therefore only practical if the time offset is known precisely (to within a fraction of the subpulse width, $\sim$a few milliseconds for B0031$-$07), or if large instantaneous bandwidths are available at a single telescope, so that the frequency evolution of the drift bands is unambiguous. Using the relative phase shifts between the three drift modes, we have deduced that if the time offset is no more than $13\,$ms (equivalent to $5^{\circ}$ in pulse longitude) from the offset obtained by aligning the average profile peaks at the two frequencies, then the emission height differences between the modes is $\lesssim 2000\,$km for mode A and $\lesssim 1000\,$km for modes B and C, or approximately a few percent of the light cylinder radius.

Within this small range of time offsets, there appears to be a greater inferred height difference generally between the drift bands of mode A than for the other two modes. If this is true, then there should also be a correspondingly greater degree of profile evolution for mode A sequences. This is indeed observed, as can be seen in both Fig. 3 and the mode-separated profiles of Smits et al. (2005). However, a more robust comparison requires a more precise determination of the emission heights than we have been able to accomplish here.

4.2.2 Subpulse separation, $P_{2}$

The fact that the drift bands bunch together at higher frequencies has led many authors to incorrectly associate the frequency behaviour of $P_{2}$ with the frequency behaviour of the average profile. However, they must be considered independently from each other owing to the fact that any frequency dependence of the former must be explained in the context of the geometry of a single field line, whereas the latter considers the (flaring) geometry of the global magnetic field structure. That the two effects are independent is evidenced by PSRs B0826$-$34 and B2020$+$28, whose $P_{2}$ appears to evolve in the opposite sense to their average profiles (Bhattacharyya et al., 2008). In both effects, however, the RFM plays an important role, since it relates the emission frequency to the local field geometry. We can expect, therefore, that a correct determination of the RFM, coupled with a correct understanding of which frequency-dependent effects dominate at different emission heights, would provide a quantitative explanation for both the single pulse and the average profile behaviour.

The connection between $P_{2}$ and the RFM is made explicit in Eq. (26), which led to the prediction of how $\overline{P}_{2}^{\prime}$ changes with frequency in Eq. (27). The functional form of (27) reveals the relationship between single pulse (i.e. $P_{2}$) and average pulse behaviour. The overall shape of the curve is set by the signs of of $d\Phi/d\varphi^{\prime}$ and $\eta$, and it is instructive to consider the possibilities. When $d\Phi/d\varphi^{\prime}<0$, the resulting curve has a vertical asymptote at a (maximum or minimum, depending on the sign of $\eta$) frequency of

[TABLE]

at which point $P_{2}$ explodes. However, no such drastic qualitative change in $P_{2}$ has ever been observed, so we can rule out this scenario. When $d\Phi/d\varphi^{\prime}>0$, Eq. (27) predicts that $P_{2}$ should change in the opposite sense to the average profile components. As noted above, this sort of behaviour has been observed (in B0826$-$34 and B2020$+$28), but it is far less common than the alternative. It is also troubling that $d\Phi/d\varphi^{\prime}>0$ implies (in the geometry of Gangadhara (2004)) that we are always seeing the subpulses on the trailing side of the fiducial point, which is unlikely to be true for most pulsars, especially those with large duty cycles. Only at the fiducial point does $d\Phi/d\varphi^{\prime}=0$, at which point the frequency dependence vanishes altogether.

Another challenge for this model is the high-frequency prediction that $P_{2}$ should continue to decrease at an ever greater rate. On this point, the erratic behaviour of $P_{2}$ for mode A as reported by Smits et al. (2007) has already been mentioned. Recent work by Ilie et al. (2019, submitted) have measured a mode B value of $P_{2}^{(\circ)}=(20^{+7}_{-5})^{\circ}$ at $\sim 1.4\,$GHz, which is notably smaller than measurements at similar frequencies shown in Fig. 7. Thus, $P_{2}$ may also behave similarly in modes A and B at high frequencies, both of which stand in stark contrast to the high-frequency prediction of the present model.

Thus, there are a number of problems with our model of how $P_{2}$ should change with frequency if AR are the dominant effects. In the particular case of B0031$-$07, the fit of Eq. (27) to the historical measurements of $P_{2}$ yields positive values for both $d\Phi/d\varphi^{\prime}$ and $\eta$. However, a positive value for the index, $\eta$, should be viewed with strong scepticism. All known models of RFM (see, e.g., Thorsett, 1992, and references therein) predict a negative index. Some pulsars are known to show unusual profile evolution, such as B1944$+$17, whose profile components move outward at higher frequencies (Kloumann & Rankin, 2010), and others whose more central components shift in phase much less (or not at all!) than the outer components. However, even in the unlikely scenario that these anomalies imply an unusual RFM, this cannot be the case for B0031$-$07. Low frequency observations of B0031$-$07 reveal that the primary profile component becomes double at $<100\,$MHz (Izvekova et al., 1993; Suleimanova & Pugachev, 2002). This implies (within the context of the carousel model) that the observed footpoints are indeed at smaller colatitudes at lower frequencies, i.e. the line of sight cuts the carousel more centrally at lower frequencies. The assumed dipolar geometry of the magnetic field requires that the emission height along inner field lines is higher than for outer field lines. We can therefore conclude that there is nothing unusual about the RFM of B0031$-$07 that requires exceptional consideration, and that B0031$-$07 is similar to most pulsars whose $P_{2}$ follows the same trend as the average profile components.

Given this, can we make sense of our measurement of a positive $\eta$? One possibility is that the RFM is a strong function of magnetic azimuth (and therefore of rotation phase), rendering Eq. (26) invalid. This idea has been explored by Thomas & Gangadhara (2007) and Thomas et al. (2010), who showed how the curvature of the field lines on the leading side is augmented by rotation effects, while on the trailing side the curvature is diminished. A more complete modelling of $P_{2}$ vs frequency for B0031$-$07 that takes rotation-induced curvature is beyond the scope of this analysis—we only note here that this effect (or some other effect that we haven’t considered) may be necessary for a complete understanding of the observed frequency-dependent behaviour of $P_{2}$ in B0031$-$07 if AR effects are the dominant cause.

Alternatively, AR effects may not be the dominant frequency-dependent effects at play here, which would naturally affect our analysis of the absolute phase shifts seen in the three drift modes (Fig. 6). If this is the case, one alternative is the effect of the sparks having finite size and shape, as studied in B0320+39 and B0809+74 by Edwards & Stappers (2003) and in B0943+10 by Bilous (2018). We do not attempt a similar analysis for B0031$-$07 here, but this could now be fruitfully attempted using the carousel configuration we have derived. In general, it would be useful to compare the spark-size model with the present AR-augmented model to see which of them dominates for a given emission geometry (i.e. carousel configuration and emission height). Alternative models, such as the fan-beam model of Dyks & Rudak (2015), may also be compared.

5 Conclusions

We have presented a new solution for the carousel configuration of B0031$-$07, showing that the subpulse behaviour of the three drift modes can be explained with a single carousel rotation rate. This is possible because of aliasing effects, which produce different apparent drift rates and $P_{3}$ values by only changing the number of sparks in the carousel. Among the solutions consistent with the drift rates (given in Table 3), we give as a representative solution a carousel rotation rate of $\overline{P}_{4}=16.4$ and number of sparks in each drift mode $[n_{A},n_{B},n_{C}]=[15,14,13]$. We have explained how first-order aliasing solutions in which the number of sparks differs between drift modes by a constant amount ($\Delta n=1$, in this case) produce $P_{3}$ values that are harmonically related, which we confirm in the case of B0031$-$07. Moreover, this carousel rotation period is much closer to the theoretical value predicted by RS than estimates that assume no aliasing.

We have extended the carousel model to include AR effects and shown how they affect the pulse phase at which subpulses appear. However, we found several problems with the model’s ability to interpret the frequency dependence of $P_{2}$, and that these problems were not limited to B0031$-$07 alone, but also to the wider pulsar population in which $P_{2}$ usually behaves in a similar way to the average profile components (smaller separation at higher frequencies). This indicates that either the AR-theory is incomplete (e.g. rotation induced curvature has not been included) or AR effects are not typically the dominant cause of frequency-dependent behaviour. In either case, the inferred emission height differences for the three drift modes ($\lesssim 2000\,$km for mode A and $\lesssim 1000\,$km for modes B and C, between $185$ and $610\,$MHz), as conservative as they are, should be nevertheless regarded with caution.

Despite these difficulties, the model describes the effect that AR must have on subpulses in the carousel model. As such, it provides a new method for estimating emission height differences that uses the frequency evolution of the drift bands, making it independent from other methods which use the average profile. If it can be shown that AR effects dominate the frequency-dependent behaviour of a given pulsar, this method provides a strong test of the carousel model, which uniquely assumes that the behaviour of the subpulses is governed by spark events located very close to the pulsar surface, thereby dictating the magnitude of phase shift that will be observed for a given emission height. Therefore, applications of this AR-augmented model to other bright pulsars which exhibit subpulse drifting and whose emission heights are known or estimated by other methods will be able to conclusively determine whether the observed phase shift is consistent with the carousel model.

The authors would like to thank A. A. Deshpande for many useful and fruitful discussions, and the anonymous referee for insightful feedback that greatly improved the paper. They acknowledge the contribution of an Australian Government Research Training Program Scholarship in supporting this research. GW thanks the University of Manchester for Visitor status. This scientific work makes use of the Murchison Radio-astronomy Observatory, operated by CSIRO. We acknowledge the Wajarri Yamatji people as the traditional owners of the Observatory site. Support for the operation of the MWA is provided by the Australian Government (NCRIS), under a contract to Curtin University administered by Astronomy Australia Limited. The GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. We acknowledge the Pawsey Supercomputing Centre, which is supported by the Western Australian and Australian Governments. Parts of this research were supported by The Centre for All-sky Astrophysics is an Australian Research Council Centre of Excellence, funded by grant CE110001020.