On the Geometry of Curvature Radiation and Implications for Subpulse Drifting

The phenomenon of subpulse drifting (Drake & Craft 1968; Backer 1973) offers unique insights into the emission physics of pulsars, if the correct interpretation can be found. One popular interpretation has its basis in the carousel model, originally proposed by Ruderman & Sutherland (1975), in which the geometry of the emitted pulsar beam is intimately connected with the geometry of a circular pattern of sparks near the stellar surface that rotate around the magnetic pole. More recently, Deshpande & Rankin (1999, 2001) developed a technique for mapping the two-dimensional emission beam that overcomes the inherent difficulty arising from the fact that a fixed observer's line of sight only ever makes a one-dimensional cut through it. This "cartographic transform" technique requires sufficiently long and stable pulse sequences⁵ in order to "fill in" the map completely enough to see the beam's global structure.

The cartographic transform involves a number of parameters, including α, the angle between the rotation and magnetic axes, ζ, the angle between the rotation axis and the line of sight, N, the number of beamlets in the carousel, and P₄, the rotation time of the carousel. If these parameters are either not known beforehand, or are coarsely estimated by other means, the cartographic transform can be used to test the viability of sets of parameters. In general, if one or more of the parameters put into the transform are slightly incorrect, the resulting polar beam pattern will be smeared and hence would be relatively devoid of structure. Then, a pulse sequence reconstructed from it will not correlate well with the original pulse sequence. However, the test is not entirely sufficient for finding the correct set of parameters as degeneracies may exist (e.g., due to aliasing of the rotating spark pattern modulo the pulsar rotation), which would have to be resolved by other means.

In the case of PSR B0809+74, Rankin et al. (2006b) found that frequent mode changes, nulls, and "memory" of the drift band phase across the nulls, all conspired to make it difficult to distinguish between a handful of solutions that all produced credible beam emission maps (see their paper, and references therein, for a full account of the difficulties encountered in their analysis). One of their solutions (α ≈ 9°, ζ ≈ 135, N = 10) produced a polarized beam map with a striking asymmetry, where the beamlets of one of the orthogonal polarizations appear to be skewed in a common azimuthal direction (see their Figure 4).

Under the assumption that the geometric parameters are correct, the authors connect this asymmetrical feature to "absorption" (Rankin et al. 2006a), which is an empirical phenomenon, currently lacking physical justification. Here, we ask whether this kind of beam asymmetry can possibly arise from purely geometric considerations once rotational effects are taken into account. In this context, the primary relevant rotational effects are aberration and retardation (hereafter, AR effects), which are both height-dependent effects that cause emission to appear at an earlier phase than otherwise expected (Blaskiewicz et al. 1991; Dyks et al. 2004b).

The treatment of aberration begins by making a few simple assumptions about how the emitted radiation would appear in the corotating frame (CF), and then using the principles of relativistic aberration to translate the Poynting vector into the observer's inertial frame (IF). Dyks et al. (2004b) show that if you assume that the Poynting vector is parallel to the tangent of the local magnetic field at the emission point, then the angle between the Poynting vector and the tangent of the local magnetic field in the IF, η_ab, is, to first order,

$$\begin{eqnarray}&&{\eta }_{\mathrm{ab}}\approx {r}^{{\prime} }\sin {\theta }_{z},\end{eqnarray} \tag{ 1 }$$

where ${r}^{{\prime} }=r/{r}_{L}$, r is the distance of the emission point from the origin placed at the stellar center, ${r}_{L}=c/{\rm{\Omega }}$ is the light cylinder radius, Ω = 2π/P is the angular frequency of pulsar rotation, and θ_z is the angle that the position vector of the emission point makes with the rotation axis $\hat{{\boldsymbol{\Omega }}}\equiv \hat{{\boldsymbol{z}}}$. The observable effect of this aberrational correction is a shift of the profile toward earlier phases,

$$\begin{eqnarray}&&{\rm{\Delta }}{\phi }_{\mathrm{ab}}\approx \displaystyle \frac{{\eta }_{\mathrm{ab}}}{\sin \zeta }.\end{eqnarray} \tag{ 2 }$$

Dyks et al. (2004b) point out, however, that the assumption of parallel vectors in the CF requires that the Lorentz factors of the emitting particles are relatively high, a requirement that is necessary in any case for curvature radiation to be in the observed range of radio frequencies.

The necessity of translating between the CF and the IF would not be required if the trajectories of the particles were already known in the IF. In this case, the Poynting vector of the curvature radiation could simply be identified with the instantaneous velocity of the particles (and of course across the width of the individual particle beam, which goes ∼1/γ). It is well established that the magnetic field lines are "frozen into" the corotating magnetosphere, and that particles are tightly constrained to move along magnetic field lines (at least, well within the light cylinder radius) as "beads on a wire" (e.g., Goldreich & Julian 1969; Sturrock 1971; Ruderman & Sutherland 1975).

The problem of determining the trajectory of a particular particle on a particular field line from first principles is a non-trivial exercise. In Thomas & Gangadhara (2007), the equations of motion for a particle in the (vacuous) magnetosphere are derived and evaluated numerically to produce particle trajectories that depend on both the functional form of the magnetic field and the magnetic inclination angle α. They found that there is significant contribution, due to rotation effects, to the curvature of field lines on the leading side of the magnetic axis as compared to those on the trailing side. Indeed, even particles on field lines very close to the magnetic axis have trajectories with significant rotation-induced curvature even though the field lines themselves have very small curvature.

Thomas & Gangadhara (2007) remind us that in the emission region, the various forces acting on a magnetospheric particle are expected to balance in such a way that the particle is tightly constrained to remain on a single field line. This is because, in the words of Radhakrishnan & Deshpande (2001), "any transverse momentum and energy would be radiated away "instantly," and the charged particles would be in their lowest Landau levels." This implies that in regions where this balance is maintained, there is a relatively simple way to deduce the velocities of particles in the emission region by geometric arguments, without needing to evaluate the competing forces acting on the particles. Indeed, a sufficient set of assumptions to derive a particle's trajectory are that (1) we know the time-dependent magnetic field in the IF (and that the observer is at rest with respect to the center of the star), (2) particles (in the emission region) move one-dimensionally along field lines as beads on a wire, and (3) we know the particle's speed in the IF, $\beta \equiv v/c$. In Section 2 we derive the particle velocities under the above assumptions, using a simple dipole field, and compare the results of this derivation to the more traditional approach via aberration as presented by Dyks et al. (2004b) and others. We then proceed in Section 3 to use the derived analytical expression of the velocity field to predict the shape of the beam, both analytically in the case of an aligned rotator, and numerically in the general case. The geometric model is then applied in Section 5 to PSRs B0809+74 and B2034+19, two subpulse drifters with contrasting drift band morphologies, to demonstrate the viability of the geometric approach. Finally, we include a discussion on the conditions required to observe various geometric effects in drift bands.

The requirement of particles to stay on a corotating field line provides a strong geometric constraint on the possible trajectories of the particle. Let ${\boldsymbol{C}}(s,t)$ be a parameterization of a given field line. Then corotation implies that ${\boldsymbol{C}}(s,t+{\rm{\Delta }}t)={R}_{z}({\rm{\Delta }}t){\boldsymbol{C}}(s,t)$, where ${R}_{z}({\rm{\Delta }}t)$ is the matrix representing constant-speed rotation about the z-axis after time Δt. An observer in the IF who measures a particle at times t and Δt will see the particle move at some non-zero angle, η, to the local magnetic field. The average speed measured in the IF is completely determined by the motion of the field line, which is known, and the angle η. Indeed, the fact that the particle's speed in the IF is necessarily v < c is equivalent to a finite range of allowed values of η, as illustrated in Figure 1. Moreover, the extreme values of η correspond to the limit $v\to c$.

Zoom In Zoom Out Reset image size

To find which values of η correspond to a given velocity, we take the vector sum illustrated in Figure 1 in the limit ${\rm{\Delta }}t\to 0$. In this limit, the vector $\overline{{{PP}}^{{\prime} }}$ becomes the particle's instantaneous velocity as measured in the IF, the vector $\overline{{PQ}}$ becomes the azimuthal velocity, and the vector $\overline{{{QP}}^{{\prime} }}$ is tangent to the local magnetic field line at P. Thus, $\overline{{{PP}}^{{\prime} }}=\overline{{PQ}}+\overline{{{QP}}^{{\prime} }}$ becomes a vector triangle for which we know two sides (the particle's measured or assumed speed and the azimuthal velocity) and an angle (between the azimuthal direction and the local magnetic field), and which therefore can be solved completely.

As implied by Figure 1, there are in general two solutions for a given particle speed. By solving the vector triangle for v ≈ c, the normalized velocity vector can be expressed as

$$\begin{eqnarray}&&\hat{{\boldsymbol{v}}}={\rho }^{{\prime} }\hat{{\boldsymbol{\phi }}}+\left[-{\rho }^{{\prime} }(\hat{{\boldsymbol{\phi }}}\cdot \hat{{\boldsymbol{B}}})\pm \sqrt{1-{\left({\rho }^{{\prime} }\right)}^{2}(1-{\left(\hat{{\boldsymbol{\phi }}}\cdot \hat{{\boldsymbol{B}}}\right)}^{2})}\right]\hat{{\boldsymbol{B}}}\end{eqnarray} \tag{ 3 }$$

where ${\rho }^{{\prime} }=\rho /{r}_{L}$ is the azimuthal corotation speed normalized to the speed of light, ρ is the point's perpendicular distance from the $\hat{{\boldsymbol{z}}}$-axis, $\hat{{\boldsymbol{\phi }}}$ is the azimuthal unit vector, and $\hat{{\boldsymbol{B}}}$ is the unit tangent to the local magnetic field in the inertial observer frame. In most cases in the region of interest (i.e., above the polar cap), the positive solution corresponds to outward-flowing particles and the negative solution to inward-flowing particles. Throughout the rest of this paper, we will consider only the positive solution. Moreover, because Equation (3) is defined at all points in the magnetosphere, it describes a unique vector field (up to the sign of the radical), from which all other geometric properties can be directly derived, as is shown later in Section 3.

The particle acceleration field can be directly computed from the velocity field. Like the magnetic field, the velocity field "corotates" with the pulsar, and in the context of a time-dependent velocity field, the acceleration field must be identified with the material derivative

$$\begin{eqnarray}&&{\boldsymbol{a}}=\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}+{v}_{x}\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial x}+{v}_{y}\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial y}+{v}_{z}\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial z}.\end{eqnarray} \tag{ 4 }$$

Because the velocity ${\boldsymbol{v}}$ has been constructed to have constant magnitude (v ≈ c), the acceleration can be reduced to

$$\begin{eqnarray}&&{\boldsymbol{a}}=c\left(\displaystyle \frac{\partial \hat{{\boldsymbol{v}}}}{\partial t}+{v}_{x}\displaystyle \frac{\partial \hat{{\boldsymbol{v}}}}{\partial x}+{v}_{y}\displaystyle \frac{\partial \hat{{\boldsymbol{v}}}}{\partial y}+{v}_{z}\displaystyle \frac{\partial \hat{{\boldsymbol{v}}}}{\partial z}\right).\end{eqnarray} \tag{ 5 }$$

Equations (3)–(5) were derived with an arbitrary magnetic field, but one must be careful to use only (time-dependent) magnetic fields that are physically meaningful in the IF. For an investigation of first-order effects, a static rotating dipole is a sufficiently good approximation, but throughout this paper, and in the accompanying numerical code,⁶ we have implemented the vacuum field derived originally by Deutsch (1955) (see Arendt & Eilek 1998; Dyks & Harding 2004).

The salient feature of the velocity field as defined by Equation (3) is that there is a unique emission direction associated with each emission location. Therefore, assuming that magnetospheric propagation effects are negligible (which may very easily not be the case, e.g., see Barnard & Arons 1986), the beam shape can be predicted for a given set of emission regions. A complete model also requires that the spectral output of each point is known, but in this analysis we will ignore any frequency dependence and only consider emission to be either on or off.

Before presenting the numerical beam shapes, we review the treatment of beam shapes in the slow rotation limit, and then derive the distortions introduced by finite rotation speeds. In this and the following sections, $(r,\theta ,\phi )$ refer to spherical coordinates aligned with the magnetic axis. That is, r, θ, and ϕ are the radial distance, the magnetic colatitude, and the magnetic azimuth, respectively.

In the absence of rotation, the velocity field is everywhere parallel to the magnetic field, as can be confirmed by letting ${r}_{L}\to \infty$ (i.e., ${\rho }^{{\prime} }\to 0$) in Equation (3). In this case, emission will be axisymmetric about the magnetic axis, $\hat{{\boldsymbol{\mu }}}$, and we can define the beam magnetic colatitude⁷ to be ${\rm{\Gamma }}={\cos }^{-1}(\hat{{\boldsymbol{B}}}\cdot \hat{{\boldsymbol{\mu }}})$. In the same limit, the Deutsch field reduces to a static dipole, and in this case, this colatitude becomes (for a fuller treatment, see Gangadhara & Gupta 2001; Gangadhara 2004)

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{nr}}={\tan }^{-1}\left(\displaystyle \frac{3\sin \theta \cos \theta }{3{\cos }^{2}\theta -1}\right),\end{eqnarray} \tag{ 6 }$$

with the well known small angle approximation Γ ≈ (3/2)θ. The subscript "nr" here indicates nonrotation. The emission, like the magnetic field, will have no azimuthal component in the magnetic coordinate system.

With rotation, the important quantity for the beam shape is the angle between the velocity field and the magnetic axis, i.e., the deflected beam colatitude ${\rm{\Gamma }}={\cos }^{-1}(\hat{{\boldsymbol{v}}}\cdot \hat{{\boldsymbol{\mu }}})$. The expression of this quantity in terms of spherical coordinates is a formidable exercise in algebra, but we can compare in more detail the more simplified case of an aligned rotator, α = 0. Then,

$$\begin{eqnarray}&&\cos {\rm{\Gamma }}=(3{\cos }^{2}\theta -1)\sqrt{\displaystyle \frac{1-{\left({r}^{{\prime} }\sin \theta \right)}^{2}}{3{\cos }^{2}\theta +1}}.\end{eqnarray} \tag{ 7 }$$

In the small angle approximation, we find

$$\begin{eqnarray}&&{\rm{\Gamma }}\approx \sqrt{{\left({r}^{{\prime} }\right)}^{2}+\displaystyle \frac{9}{4}}\,\theta +O({\theta }^{3}),\end{eqnarray} \tag{ 8 }$$

showing that the beam is dilated by an extra factor that depends on the emission height.

In the context of a given field line, r' and θ are coupled. For a dipole geometry, the relationship between them is: $r\propto {\sin }^{2}\theta$, where the constant of proportionality is the field line's maximal distance from the origin. Here, we follow the convention of Gangadhara & Gupta (2001) and others, and identify a field line by the magnetic colatitude of the point at which it pierces the pulsar surface, i.e., the "footpoint" colatitude θ_p. Moreover, we normalize the footpoint colatitude to the colatitude of the last open field lines, ${\theta }_{L}={\sin }^{-1}\sqrt{{r}_{p}/{r}_{L}}$, so that the quantity $s\equiv {\theta }_{p}/{\theta }_{L}$ equals 0 at the magnetic pole and unity at the polar cap radius. By fixing s, we can make the substitution

$$\begin{eqnarray}&&{r}^{{\prime} }\to \displaystyle \frac{{r}_{p}{\sin }^{2}\theta }{{r}_{L}{\sin }^{2}{\left(s{\theta }_{L}\right)}^{2}},\end{eqnarray} \tag{ 9 }$$

in which case the approximation (8) reduces to Γ ≈ (3/2)θ as before, with the s term only making an appearance in the 5th order term in θ.

As well as a slight expansion, the beam also acquires an azimuthal component, which due to axisymmetry also depends only on r' and θ. If the intensity of the pulsar beam is also axisymmetric, the azimuthal component becomes irrelevant (for the aligned rotator), since it only serves to rotate the beam around the common axis. However, the intensity profile generated by a rotating carousel is (or may be) only axisymmetric after averaging over a carousel rotation. In order to appreciate the effects of aberration on the global beam shape, we momentarily neglect the effects of carousel rotation by considering the slow rotation limit (i.e., ${P}_{4}\to \infty$). In this limit, each beamlet originates from a narrow tube of field lines emerging from a discrete spark at the tube's base. The beam shape due to a single narrow tube will in general be some form of spiral, as emission from different heights along the tube is deflected by different amounts (Maan & Deshpande 2014). Of course, the whole spiral will not necessarily be present because the emission is presumably only generated within a finite range of heights at a given frequency. Nevertheless, for a given line of sight cut, it may be the case that a large enough range of heights are sampled that the map produced by the cartographic transform reveals the spiral nature of the individual beamlets.

To find the beam's azimuthal component (in the aligned case), we first compute the quantity ${\rm{\Delta }}{\phi }_{\mathrm{ab}}={\tan }^{-1}(-{v}_{y}/{v}_{x})$. Axisymmetry allows us to calculate the azimuthal aberration by simply evaluating Δϕ_ab in the xz-plane (i.e., when ϕ = 0). This gives

$$\begin{eqnarray}&&\tan ({\rm{\Delta }}{\phi }_{\mathrm{ab}})=-\displaystyle \frac{{r}^{{\prime} }}{3\cos \theta }\sqrt{\displaystyle \frac{3{\cos }^{2}\theta +1}{1-{\left({r}^{{\prime} }\sin \theta \right)}^{2}}}.\end{eqnarray} \tag{ 10 }$$

The above result may be more simply derived by noting that in Figure 2, the common sides of the triangles formed by ${\boldsymbol{v}}$, ${{\boldsymbol{v}}}_{\phi }$, and ${{\boldsymbol{v}}}_{B}$ must be related by

$$\begin{eqnarray}&&\cos ({\rm{\Delta }}{\phi }_{\mathrm{ab}})=\displaystyle \frac{\tan {{\rm{\Gamma }}}_{\mathrm{nr}}}{\tan {\rm{\Gamma }}},\end{eqnarray} \tag{ 11 }$$

into which Equations (6) and (7) can be substituted.

Zoom In Zoom Out Reset image size

Having obtained a complete description of the deflected emission beam, we must now incorporate the effects of retardation into the analysis. Similarly to Γ above, we can define the retardation angle, Δϕ_ret, as a valid quantity at every point in the magnetosphere, regardless of whether or not the observer is in a direction to observe it:

$$\begin{eqnarray}&&{\rm{\Delta }}{\phi }_{\mathrm{ret}}=-{r}^{{\prime} }(\hat{{\boldsymbol{r}}}\cdot \hat{{\boldsymbol{v}}}),\end{eqnarray} \tag{ 12 }$$

which, for an aligned rotator becomes

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{\phi }_{\mathrm{ret}} & = & -2{r}^{{\prime} }\cos \theta \sqrt{\displaystyle \frac{1-{\left({r}^{{\prime} }\sin \theta \right)}^{2}}{3{\cos }^{2}\theta +1}}\\ & = & -\displaystyle \frac{2{r}^{{\prime} }\cos \theta \cos {\rm{\Gamma }}}{3{\cos }^{2}\theta -1}.\end{array}\end{eqnarray} \tag{ 13 }$$

The negative sign serves the same purpose as that in Equation (10), namely, to shift the observed emission to earlier rotation phases. This shift is equivalent to a further azimuthal distortion of the overall beam pattern, which in the aligned case serves to enhance the curvature of the spiral shape of the beam. For a dipolar field, Equation (12) reduces to the approximation ${\rm{\Delta }}{\phi }_{\mathrm{ret}}\approx -{r}^{{\prime} }+O({{r}^{{\prime} }}^{3})$, in agreement with Dyks et al. (2004a).

By choosing a fixed s, the equations above constitute a complete description of the spiral beam pattern, as illustrated in Figure 3. The foregoing construction is also valid (albeit algebraically daunting) for pulsars of arbitrary inclination angle α.

Zoom In Zoom Out Reset image size

Emission from field lines with smaller s suffer from greater azimuthal distortion. This makes intuitive sense: the curvature of these field lines is smaller, and so the particles must climb to greater heights in order to produce beams of equivalent apparent magnetic colatitude. In the case of any one particular pulsar, however, it generally remains unknown exactly which magnetic field lines are the active ones. On the other hand, if azimuthal distortion of the beam can be measured, and assuming there are no significant additional distorting effects, the footpoints colatitudes can be constrained.

A true comparison of the predictions of the present model and the polar maps generated by the cartographic transform requires the simulation of pulsestacks, which can then be transformed. In this section, we describe in detail the process of producing simulated pulsestacks from the geometric model outlined above.

Clearly, the geometric framework alone is insufficient to simulate observed pulsestacks; it can only say when emission from a given location in the magnetosphere would be seen by the observer, not how bright that emission would be at a given frequency. However, the primary consideration here is the morphology of the drift bands, i.e., when the subpulses in successive pulses arrive in relation to each other. By virtue of the underlying assumption of the carousel model—that the emissivity of a field line is directly related to the spark activity at its base—we can use a simple model for estimating the observed intensities from specific emission locations in which they are proportional to the spark activity at their respective footpoints.

For an observer whose line of sight makes an angle ζ to the rotation axis, the only points on a given magnetic field line that are potentially visible are those for which

$$\begin{eqnarray}&&\hat{{\boldsymbol{v}}}\cdot \hat{{\boldsymbol{\Omega }}}=\cos \zeta \end{eqnarray} \tag{ 14 }$$

(in the aligned case, this condition is equivalent to Γ = ζ). When ζ > α, there is a unique height on each magnetic field line above the polar cap where this condition is satisfied. On the other hand, when ζ < α, field lines on the equatorial side of the magnetic axis start at the surface with associated velocities that make an angle to the rotation axis >α, and this angle only increases at greater emission heights as the field line drops away, never becoming equal to ζ. Then $\zeta \lt \alpha \lt {\cos }^{-1}(\hat{{\boldsymbol{v}}}\cdot \hat{{\boldsymbol{\Omega }}})$, so Equation (14) can never be satisfied and these field lines never come into view of the observer. Field lines on the poleward side, however, initially curve "upwards" and then "backwards," so that the angle ${\cos }^{-1}(\hat{{\boldsymbol{v}}}\cdot \hat{{\boldsymbol{\Omega }}})$ starts at some value <α, decreases through ζ until some minimum value is reached, and then increases again as the magnetic field drops away on the other side. These field lines therefore include visible points at two different heights, whose exact values can be found by solving Equation (14) numerically. Since the second visible point is typically much higher than the first and appears near the anti-fiducial point, they are not expected to be observed in type S_d pulsars (i.e., whose profiles consist of a single conal component) and we hereafter consider only the lower visible point.

The fact that each visible magnetic field line is only observed at most at one height implies that any arbitrary intensity profile or pulsestack can be obtained by a careful choice of (dynamic) spark pattern on the surface. However, we note that multiple field lines can contribute to the emission observed at a given phase, so that at all phases a range of heights is (potentially) observed. In this geometric model, having a range of emission heights is necessary to produce a spread in the polarization angle at a given rotational phase, because just as each point in the magnetosphere has associated with it a unique velocity from Equation (3), each point also has a unique acceleration vector from Equation (4), and therefore a unique polarization angle. Thus, a realistic model requires that a full, two-dimensional spark pattern is used, in order that a range of heights and therefore a range of polarization angles are sampled at each phase. However, for ease of computation, the simulated pulsestacks presented here are generated from a one-dimensional carousel at a fixed value of s whose "sparks" have Gaussian profiles in magnetic azimuth. By choosing only a single s value and making the rest of the polar cap "inactive," we are forcing the polarization angle to adopt a unique value at each phase—that is, every pulse will have the exact same sweep of polarization angles—which we will interpret as the "average" polarization angle.

The steps to create a pulsestack are then as follows. For each footpoint in the (one-dimensional) carousel, the field line is traversed numerically using Runge–Kutta integration until a point is found that has the property that $\hat{{\boldsymbol{v}}}\cdot \hat{{\boldsymbol{\Omega }}}=\cos \zeta$, where $\hat{{\boldsymbol{v}}}$ is calculated from Equation (3). This is the so-called "visible point" of that field line. The polarization angle associated with the visible point is then calculated by projecting the acceleration vector of Equation (4) onto the observer's sky plane (i.e., normal to $\hat{{\boldsymbol{v}}}$). This is done to verify that the chosen parameters are consistent with the observed (average) polarization swing.

The observed phase is determined by calculating Δϕ_ab + Δϕ_ret. The first term gives the rotation phase at which the emission must have taken place, and the second term is the retardation shift due to photon flight times. Next, we find the rotation phase at which a spark event at the field line's footpoint would have occurred in order that the resulting particle stream (assumed to be traveling at ultra-relativistic speeds) would reach the visible point at phase Δϕ_ab,

$$\begin{eqnarray}&&{\rm{\Delta }}{\phi }_{\mathrm{sp}}={\rm{\Delta }}{\phi }_{\mathrm{ab}}-{{\ell }}^{{\prime} },\end{eqnarray} \tag{ 15 }$$

where ${{\ell }}^{{\prime} }\equiv {\ell }/{r}_{L}$, and ℓ is the distance traversed by a particle as it corotates with the pulsar and climbs from the surface to the emission point, which is calculated numerically at the same time that the visible point is being found. For pulse number p in the pulsestack, Δϕ_sp is converted (via the pulsar rotation period, P₁) into a time ${t}_{\mathrm{sp}}={P}_{1}(p+{\rm{\Delta }}{\phi }_{\mathrm{sp}}/2\pi )$, which can then in turn be converted into an intensity via the carousel model

$$\begin{eqnarray}&&I({t}_{\mathrm{sp}},\phi )\propto \displaystyle \sum _{n=1}^{N}\exp \left[-\displaystyle \frac{{\left(\phi -2\pi \left(\tfrac{n}{N}+\tfrac{{t}_{\mathrm{sp}}}{{P}_{4}}\right)\right)}^{2}}{2{\sigma }^{2}}\right],\end{eqnarray} \tag{ 16 }$$

where ϕ is the magnetic azimuth of the footpoint and σ parametrizes the angular width of the individual spark profiles.

4.1. Phase-dependent Intensity Modulation

The foregoing procedure will produce a pulsestack where emission is "observed" across all 360° of longitude, regardless of the local conditions at the emission point, as illustrated in Figure 4. In other words, it follows the assumption of the carousel model that sparks activity at the footpoint of a field line is a necessary, but not a sufficient condition for coherence to occur on that field line. Some extra condition, presumably connected to the emission mechanism, must therefore be imposed in order to decide how the drift bands are modulated over the pulse phase.

Zoom In Zoom Out Reset image size

If a mechanism such as particle bunching is assumed (as first suggested by Ruderman & Sutherland 1975, but see, e.g., Melrose 2017 for arguments against bunching as a viable mechanism), then the beam intensity is proportional to N², where N is the number of particles in a "bunch" (i.e., within a volume element with its linear dimension smaller than the emitted wavelength in the co-moving frame). In this case, the observed emission pattern for curvature radiation from a single particle can be calculated, and scaled up by the square of the local plasma density at each emission site. Thomas et al. (2010) describe one possible method for doing such a calculation, although it should be noted that they incorporate the phase-shifting effects of aberration and retardation post hoc. To our knowledge, no similar method exists for estimating intensity profiles assuming other possible mechanisms, such as maser mechanisms (Ginzburg & Zhelezniakov 1975) or plasma wave instabilities (Cheng & Ruderman 1977).

It should be realized, however, that the positions of the drift bands in observed rotation phase are, under the assumptions of the carousel model, agnostic to which emission mechanism is ultimately adopted. If, therefore, the goal is to study how the drift bands are affected by the geometry of the magnetosphere, one can avoid (or at least defer) the unenviable task of choosing and modeling one particular coherent emission mechanism at the expense of the others. If a comparison with real data is still needed, it suffices to modulate the simulated pulsestack with an actual average profile at the desired frequency, which would include asymmetries, if any, resulting from both the emission physics and geometrical considerations.

5.1. Drift Band Morphology

5.2. Appearance of Reconstructed Beamlets

The simulated pulsestacks can be subjected to all the standard analyses used to study subpulse drifting. In particular, they are amenable to the reconstruction of polar beam maps via the cartographic transform of Deshpande & Rankin (2001). As a case study, we investigate one of the viewing geometries proposed by Rankin et al. (2006b) for PSR B0809+74, whose transformed beam map shows a remarkable azimuthal skewing of the beam maps (see their Figure 4). The viewing geometry assumed α = 9°, ζ = 135, and number of sparks N = 10. Even though the authors emphasize the difficulty of distinguishing between many sets of viable parameters, they argue that in this case aliasing is not present, in which case P₃ ≈ 11.1 P₁ is the time interval between successive drift bands at a fixed rotation phase, giving ${P}_{4}=N\times {P}_{3}\approx 143\,{\rm{s}}.$

One parameter remains free: the (normalized) radius of the carousel, s. We first tested which values in the range 0 < s < 1 produced polarization angle curves consistent with the observations (for one of the orthogonal polarization modes only). Figure 5 shows the polarization angle curves predicted by the geometric model. As expected, values for relatively large s agree better, as these correspond to relatively low emission heights, at which the geometric model more closely approximates the rotating vector model (RVM), which was employed in the determination of the viewing geometry parameters. In particular, we note that values in the range 0.4 ≲ s ≲ 1.0 agree sufficiently well with the polarization angles in the pulse window.

Zoom In Zoom Out Reset image size

After choosing a value for s, the pulsestack can be simulated according to the procedure described in Section 4. Figure 6 shows a simulated intensity pulsestack for pulsars B0809+74 and B2034+19, using parameters reported in Rankin et al. (2006b) and Rankin (2017), respectively. The average profile used to modulate the pulsestack was a ∼150 MHz profile retrieved from the European Pulsar Network database,⁸ supplied by Noutsos et al. (2015). The authors assigned the fiducial point in their profile by means of an RVM (Radhakrishnan & Cooke 1969; Komesaroff 1970) fit to the observed polarization curve. Since (by virtue of Figure 5) the simulated polarization curve is sufficiently close to the RVM for the range 0.4 < s < 1.0, the association with their empirically determined fiducial point with the simulated fiducial point is justified. Varying s within this range does not significantly alter the shape of the drift bands.

Zoom In Zoom Out Reset image size

Having obtained a simulated pulsestack, we then applied the cartographic transform in order to obtain a simulated polar map of the pulsar beam. Figure 7 shows the results for two choices of footprint polar radius, s = 0.5 and s = 0.25. When s values are chosen in the range 0.4 ≤ s ≤ 1.0, no significant azimuthal distortions were observed, but smaller values of s yield very noticeable azimuthal distortions with similar characteristics as that observed by Rankin et al. (2006b). Curiously, the direction of the azimuthal distortions differ, with the simulated beamlets being skewed in the counter-clockwise (CCW) direction, while the real data is skewed in the clockwise (CW) direction. This is discussed further in the next section.

Zoom In Zoom Out Reset image size

6.1. Comparison of Simulated and Real Data

6.2. Effect on the Interval between Successive Subpulses (P₂)

An important feature of the geometric model is the inclusion of the time taken for the spark information to reach the visible point, a non-negligible time, which to our knowledge is entirely lacking in other models of carousel behavior. The quantity P₂, which measures the elapsed time between the observation of two adjacent sparks, must therefore include not only AR effects and continuous carousel rotation (affecting the apparent spacing between the sparks, a theme discussed in Yuen et al. 2016), but also the difference of path lengths of the respective particle trajectories. The length of the trajectories traced by particles in the IF (ℓ') are necessarily different from the path length along the magnetic field lines on which they reside—rotation effects should not be neglected.

As described above, the trajectory lengths are calculated numerically, but we can get a sense of the magnitude and behavior of the effect by taking the first-order approximation ℓ' ≈ r'. Then the interval between the arrival of successive subpulses will depend on how far the carousel rotates in time ∼Δr/c, which in radians is $\tfrac{{\rm{\Delta }}{r}^{{\prime} }{P}_{1}}{{P}_{4}}$. On the leading side of the fiducial point, the emission heights start large, approach a minimum near the fiducial point, and grow large again on the trailing side (see Figure 8). The inclusion of trajectory lengths therefore either serves to artificially augment or diminish the observed value of P₂ during leading phases (depending on whether the line of sight cuts across the beam pattern in the same or opposite directional sense, respectively, as the carousel rotation) and to have the opposite effect on P₂ on the trailing side. This is similar to retardation, except that retardation always compresses the observed emission on the trailing side (Dyks et al. 2010). Therefore, the present effect can either enhance or suppress the effect of retardation, depending on the relative signs of P₁, P₄, and β = ζ − α.

Zoom In Zoom Out Reset image size

6.3. Effect of Finite Spark Size

6.4. The Rotational Asymmetry of B0809+74's Reconstructed Beam

6.5. Other Considerations

We have described a geometric model for pulsar emission that assumes coherent curvature radiation from a corotating, relativistic plasma, consistent with AR effects. The further assumption (implicit in the carousel model) that the "activity" along the field lines is only dependent on the spark activity at their footpoints near the stellar surface enables the model to predict the morphology of drift bands without reference to any particular coherent emission mechanism.

We have explored two specific pulsars, PSRs B0809+74 and B2034+19, with contrasting viewing geometries (nearly-aligned and oblique rotators, respectively) and shown that the geometric model can reproduce the essential characteristics of the observed pulsestacks, including the modulation periodicity P₂, and drift band curvature. In the case of B0809+74, we have shown how the geometric model is capable of generating rotational asymmetries in the polar maps produced by the cartographic transform, qualitatively similar to those observed by Rankin et al. (2006b). However, in this particular case, the directional sense of the asymmetries and the emission height range needed to reproduce them in simulation are difficult to reconcile with their interpretation in the current geometric context.

The geometric model presented in this paper is, we believe, a necessary consequence of the assumptions of curvature radiation, and is therefore applicable to all radio pulsars for which some form of coherent curvature radiation is believed to be the primary emission mechanism, including but not limited to the class of subpulse drifters. The most telling tests of the validity of this model require high time resolution, broadband observations of individual pulses, especially where microstructure can be resolved—a kind of observation that is scarce, due to the difficulty of coordinating multiple telescopes for a wide spectral coverage. However, with several broadband instruments coming online (e.g., the ultra-wideband receiver at Parkes, the upgraded GMRT, and the RRI-GBT multiband receiver; Maan et al. 2013), in-depth investigations of these subtle effects will become possible in the near future.

The derivation of Equation (8) and, in general, any Taylor expansion about the magnetic axis that depends on Equations (3) and (4), requires the conversion between the observer coordination system and the magnetic coordinate system.

The observer coordinate system is defined as follows. The origin is placed at the center of the pulsar. The pulsar's rotation axis is identified with the z-axis. The x-axis is chosen to lie in the plane spanned by the z-axis and the line of sight, $\hat{{\boldsymbol{n}}}$, such that positive x points toward the distant observer (i.e., $\hat{{\boldsymbol{n}}}\cdot \hat{{\boldsymbol{x}}}\geqslant 0$). Finally, $\hat{{\boldsymbol{y}}}=\hat{{\boldsymbol{z}}}\times \hat{{\boldsymbol{x}}}$.

The magnetic coordinate system is related to the observer coordinate system by a rotation about the y-axis by the magnetic inclination angle α, followed by a rotation about the z-axis by the rotation phase angle ϕ. This brings the z-axis into alignment with the magnetic axis, $\hat{{\boldsymbol{\mu }}}$. When ϕ = 0, the magnetic axis is in the plane spanned by the x- and z-axes. We note that although the magnetic coordinate system does not represent an IF, the magnetic field (${\boldsymbol{B}}$), the velocity field ${\boldsymbol{V}}$, and the acceleration field (${\boldsymbol{A}}$) are all static in the CF. That is to say, an inertial observer would find that these fields all corotate with the pulsar, so that a measurement of the entire fields at two different times would be related by a simple rotation about the z-axis (see also Appendix A of Dyks & Harding 2004).

Therefore, there is no loss of generality by evaluating these fields at a single rotation phase, say at the fiducial point (ϕ = 0) and rotating it a posteriori, or equivalently letting the line of sight "rotate" in the opposite direction by the same amount. In many cases, this makes the evaluation of Equation (3) and other derived quantities much simpler since they can be expressed independently of ϕ. Furthermore, because the only radiation that would be observed is that produced by particles whose velocity vector is parallel to the line of sight, we can simply replace all instances of $\hat{{\boldsymbol{n}}}$ with $\hat{{\boldsymbol{v}}}$ and tacitly ignore times and locations in the magnetosphere where they are not parallel.

The rotation matrices are given many times in the literature, and are not repeated here. For example, the Appendix of Yuen et al. (2016) includes the relevant transformation matrices, as well as the matrix used to convert between Cartesian and spherical polar coordinates, which is also used in this work for Taylor expansion about the magnetic axis (i.e., about θ = 0).

This research was conducted in part by the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and in part by an Australian Government Research Training Program Scholarship. G.W. thanks the University of Manchester for Visitor status. The authors would also like to thank the anonymous referee for some insightful discussion and suggestions, which improved the presentation and clarity of the paper.