A&A 393, 733-748 (2002)
DOI: 10.1051/0004-6361:20021067
R. T. Edwards1 - B. W. Stappers1,2
1 - Astronomical Institute ``Anton Pannekoek'',
University of Amsterdam,
Kruislaan 403, 1098 SJ Amsterdam, The Netherlands
2 -
Stichting ASTRON, Postbus 2, 7990 AA Dwingeloo, The Netherlands
Received 27 November 2001 / Accepted 19 July 2002
Abstract
The basic form of drifting sub-pulses is that of a periodicity whose
phase depends (approximately linearly) on both pulse longitude and
pulse number. As such, we argue that the two-dimensional Fourier
transform of the longitude-time data (called the Two-Dimensional
Fluctuation Spectrum; 2DFS) presents an ideal basis for studies of
this phenomenon. We examine the 2DFS of a pulsar signal synthesized
using the parameters of an empirical model for sub-pulse behaviour.
We show that the transform concentrates the modulation power to a
relatively small area of phase space in the region corresponding to
the characteristic frequency of sub-pulses in longitude and pulse
number. This property enables the detection of the presence and
parameters of drifting sub-pulses with great sensitivity even in data
where the noise level far exceeds the instantaneous flux density of
individual pulses. The amplitude of drifting sub-pulses is modulated
in time by scintillation and pulse nulling and in longitude by the
rotating viewing geometry (with an envelope similar to that of the
mean pulse profile). In addition, sub-pulse phase as a function of
longitude and pulse number can differ from that of a sinusoid due to
variations in the drift rate (often associated with nulling) and
through the varying rate of traverse of magnetic azimuth afforded by
the sight line. These deviations from uniform sub-pulse drift
manifest in the 2DFS as broadening of the otherwise delta-function
response of a uniform sinusoid. We show how these phase and amplitude
variations can be extracted from the complex spectrum.
Key words: methods: analytical - pulsars: general
Very soon after the discovery of pulsars, it was noticed that the intensity data sometimes showed "second periodicities'' within each pulse, and that the relative phase of this periodic signal drifted in an organised manner from pulse to pulse (Drake & Craft 1968). The striking patterns made by drifting sub-pulses in two-dimensional pulse longitude-time diagrams seemed to be saying something about the fundamentals of radio pulsar emission, but the question was, and still is, "what?'' Numerous studies of the phenomenon were conducted in the first few years of pulsar science (e.g. Taylor et al. 1969; Sutton et al. 1970; Cole 1970; Backer 1970a,b; Taylor & Huguenin 1971; Backer 1973; Page 1973), their frequency steadily decreasing as the limitations of existing online (recording) and offline (analysis) equipment were reached. Along with attempts at fitting the data with sub-pulses and tracking their drift, the Fourier technique developed by Backer in the references above (now known as the longitude-resolved fluctuation spectrum, or LRFS) came to be perhaps the most widely used means of examining the average properties of sub-pulse drift.
With the availability of modern digital back-ends and convenient offline computing power, we believe that there is now potential for renewed progress in the field of drifting sub-pulses. We describe in this paper a means of analysis of drifting sub-pulse data using what may be considered an extension of the Fourier work of Backer. The technique was first conceived as a means of detecting and characterising drifting sub-pulses in the large population of moderately weak pulsars that have been discovered since the 1970s when most studies took place. It is also useful, when sufficient signal-to-noise ratio is available, for the investigation of what might be considered the fine details of the phenomenon: deviations from constant amplitude sub-pulses with purely uniform phase evolution. We begin with a mathematical description of a drifting sub-pulse signal as a function of pulse longitude and pulse number (Sect. 2) and examine its two-dimensional Fourier response (Sect. 3). We then describe a technique for measuring the deviations of a given signal from a uniform pure sinusoid (Sect. 4.1), examine what is expected under the usual model of drifting sub-pulses (Sect. 4.2), and use this model as a basis for simulations testing the applicability of the methods of analysis described here (Sect. 4.3).
The characteristic form of drifting sub-pulses is difficult
to miss in a suitable plot of high-quality data. Plotted as intensity
as a function of pulse longitude and pulse number, a distinctive
"tire tread'' pattern is made as the sub pulses drift monotonically
in longitude with each successive pulse. The form is essentially that
of a windowed two-dimensional sinusoid. The characteristic period of
the sub-pulses in pulse longitude (P2) and pulse number (when expressed as a time interval
) may be measured by a
variety of techniques and is usually the first step in any analysis of
pulsars exhibiting drifting sub-pulses.
Naturally the actual signal is not a uniform, pure two-dimensional sinusoid. The intensity at a given instant depends on the pulse longitude (the overall emission is, after all pulsed at the spin period P1) and also has a time variability that may be intrinsic to the pulsar or induced by interstellar scintillation. The spacing of the sub-pulses may also vary in a systematic manner with time and/or pulse longitude. We shall show later (Sect. 4.2) that under the most popular model for the origin of drifting sub-pulses, we may to a very good approximation model these variations in sub-pulse spacing as a dependence of P2 on pulse longitude and/or a time-variability of P3. Observational results showing that (in specific pulsars) P2 does not vary with time (Page 1973; Unwin et al. 1978; Lyne & Ashworth 1983) and that P3 is independent of pulse longitude (Backer 1970a,b; Krishnamohan 1980; Wright 1981; Davies et al. 1984; Biggs et al. 1987) serve to confirm this assertion. The only reported exception, to our knowledge, is the sub-pulse modulation reported in the putative pulsar remnant of Supernova 1987A (Middleditch et al. 2000), which, unlike other pulsars with drifting sub-pulses, is probably caused by precession. We therefore proceed to describe the signal mathematically on this basis.
The most obvious difference between a true sub-pulse signal and a
purely periodic signal its longitude-dependent amplitude windowing. It
is expected that the sub-pulses will be scaled in amplitude as a
function of pulse phase with a shape similar to the average pulse
profile. We denote this function
,
where
is
the pulse longitude.
We may also expect some dependence of longitudinal sub-pulse
spacing on pulse longitude. Expressed another way, the sub-pulses of a
given pulse are periodic in a non-linear function of longitude. The
phase must be at least monontonic and most likely almost linear, or we
would not consider the sub-pulses to be (almost) periodic. We may
therefore treat "sub-pulse phase'' ()
in a given pulse as the
sum of a linear term and a longitude-dependent "phase envelope''
(
).
The sub-pulse signal for a given pulse can then be represented as
![]() |
(1) |
Since the deviation from pure periodicity is represented by a
longitude-dependent
amplitude scaling and phase shift, we may simply define a complex
modulation envelope
,
writing the varying sub-pulse
component of a given pulse as
![]() |
(2) |
In addition to the longitude dependence covered in the
previous section, drifting sub-pulse signals vary on time-scales
longer than one pulse period (P1). We define a time-dependent
amplitude scaling
which gives rise to the
pulse-specific A' factor used earlier. To describe the
time-dependence of the spacing of sub-pulses in time, we may take a
similar approach and define
as the difference in phase
between the true signal and a uniform sinusoid of period
.
This function is the source of the pulse-specific
terms of
the previous section. Again we can define a complex envelope
.
From the preceding two sections we are able to assemble a model for
the observed intensity of pulsar signals as a function of time and
pulse longitude. Under the assumption that the
time-dependent and longitude-dependent modulations are independent,
we may model the overall modulation as the product of the derived
longitude and time windowing:
.
Hence the model of the sub-pulse component is
The two-dimensional Fourier transform of a function f(x, y) is
defined as follows:
![]() |
(8) |
For a discretely, uniformly sampled function with
samples,
the two-dimensional
Discrete Fourier Transform (DFT) may be written as
F(u, v) | = | ![]() |
(9) |
When the signal of a pulsar with drifting sub-pulses is displayed as a
function of pulse number and longitude, the patterns made by the
sub-pulses and their drifting are very similar to those of the inverse
Fourier transform of a delta function. This fact alone suggests that
computing the (discrete) Fourier transform of the two dimensional signal
may be an interesting exercise.
Since the signal is real-valued, the result of the Fourier transform
(
)
will obey the symmetry
,
where the superscript * denotes
the complex conjugate.
From Eqs. (4), (5) and (6) the result of this operation (which we call the
Two-Dimensional Fluctuation Spectrum or 2DFS) is
![]() |
= | ![]() |
|
= | ![]() |
(10) |
![]() |
= | ![]() |
(11) |
= | ![]() |
(12) |
So we see that the 2DFS consists of three main components: one for the
response of the steady part of the pulsar emission, and two
mirror-image components for the sub-pulse modulation. The steady
component is centered at DC (
), whilst the
sub-pulse components are centered approximately at
(cycles per radian of pulse longitude),
(cycles
per time unit) where
.
Since the modulation envelopes
are not purely real-valued, their Fourier transforms are not expected
to be entirely symmetrical, and there is no unique "center'' of the
components, however in most cases the phase rotation of the envelope
will be quite minor and a peak around the position specified will be
seen.
The overall width and shape of the components depends on
the nature of the modulation envelopes. For a pulsar with no nulling
or scintillation, the steady component will have zero width
in the
direction
, and a one-dimensional inverse Fourier transform of
in the
dimension yields the average pulse
profile. If nulling and/or scintillation are significant, there may be
some broadening in the
direction however very
frequent nulling or scintillation on very short time scales would be
required to produce a strong, significantly extended component.
In addition to the effects of nulling and scintillation,
the sub-pulse components are broadened due to phase (or P3)
variation over the course of the observation.
Broadening of the components in the
axis derives
from the finite width of the amplitude envelopes (
and
)
and, in the case of the sub-pulse components, due to
phase (or P2) variation across the pulse.
The width of this component is
generally quite large due to the small duty cycle of pulsars, and will
tend to make the spectrum significantly non-zero at
,
if the number of sub-pulses observable in each
pulse is less than two.
As noted earlier, the frequency extent of the DFT is limited by the inverse of the sampling interval(s). Co-efficients for frequencies outside the range may be computed but will always be identical to those at some (calculable) point in the range due to aliasing. In the case of the 2DFS, the longitude sampling may be chosen arbitrarily (within the limits of the hardware and the bandwidth of the received radiation), but the "time'' or "pulse number'' sampling is naturally set at a value of P1. For this reason, drifting sub-pulse components outside the range -0.5 < P1/P3 < 0.5 are indistinguishable from those within it.
It is readily apparent from Eq. (6) that the
discrete two-dimensional Fourier transform can be computed by first
computing a one-dimensional DFT along each column of the input data,
and then computing the one-dimensional DFT along each row of the
result (or alternatively, along rows first and columns second). The
result of the first half of this operation when performed on the
pulsar signal
is known as the longitude-resolved
fluctuation spectrum (LRFS; Backer 1970a,b).
The expected form of this spectrum is easy to understand in the light
of the preceding discussion. Each column (i.e. line of constant
)
in the spectrum will have three components, as in the
2DFS. The steady component simply appears as a DC term equal to
,
whilst the two components for the sub-pulse
modulation appear at
.
Calculating the spectrum
of both components (denoted I1 and I2), we find:
![]() |
(13) |
![]() |
(14) |
Due to the presence of a component in both the positive and negative
half of the spectrum, it is impossible to determine the sign of
from the LRF power spectrum. Signals are effectively
aliased to the range
.
Only by examining
the sense of phase rotation of the complex co-efficients as a function
of longitude and comparing this to the sign in
of the component
being considered can the sense of the drift be determined.
Since the phase relation between longitude bins is expected to be quasi-periodic, it makes sense to perform a Fourier transform across each row of the complex spectrum in order to determine the sense of the drift and its period (i.e. P2/P1), and to concentrate the sub-pulse power in a smaller region of phase space for greater signal-to-noise ratio in P3. The result of this operation is the two-dimensional DFT of the signal, in other words the 2DFS.
Recently a new technique for sub-pulse analysis has been devised
by Deshpande & Rankin (2001). It involves the computation of the so-called
Harmonic-Resolved Fluctuation Spectrum (HRFS). Working from the
original time-series i(t), they compute the one-dimensional spectrum
and "stack'' it about the
harmonics of the signal. Thus,
![]() |
(15) |
We note that the 2DFS is in fact intimately related (after scaling and rotation) to the HRFS. The former is formed by stacking the time-domain data about the pulse period P1 and then Fourier transforming, whilst the latter is formed by Fourier transforming and stacking the result about 1/P1. That the two should be essentially identical is not immediately obvious but is easy to show (see Appendix A). In fact, the decomposition of the one-dimensional DFT into a two-dimensional DFT is the basis of parallel FFT codes (Cooley & Tukey 1965).
We note that values of
giving rise to
are
aliased to
in the 2DFS (i.e. the opposite drift sense),
however this is simply a matter of convention: by the same notion,
signals which truly show the opposite drift sense (that is, negative
P3P2) are aliased in the HRFS into the region
.
In our view, given that there is no reason to prefer one sense
of intrinsic sub-pulse drift over another, specification of
in the interval
[-0.5,0.5] makes more sense than taking
a preferred drift direction and mapping all drift to this interval
(e.g. [0:1]).
Whilst the two spectra are equivalent, a knowledge of their basis in the two-dimensional Fourier transform of the longitude-time data is of value in further analysis, as we shall see below.
With an awareness of the mathematical basis of the 2DFS in relation to the form of the drifting sub-pulse signal described in the preceding sections, a range of analysis strategies can be developed for the investigation of particular aspects of the sub-pulses.
The technique was originally conceived with the aim of detecting the presence and basic parameters (P2 and P3) of sub-pulse drift in pulsars with limited instantaneous signal-to-noise ratio. The sensitivity of this method depends only on the coherence of the drift and the ratio of its amplitude to that of any non-drifting signal. The significance of a detection can by improved arbitrarily (for stable drifters) by simply integrating more pulses, unlike in traditional single-pulse studies where each individual sub-pulse must be detectable above the noise. This makes it well-suited to studying the sub-pulse properties of a large sample of weaker pulsars.
The complex 2DFS is also of use in studying the details of sub-pulse emission in bright pulsars. In this section we show how the complex modulation envelopes can be extracted from the spectrum, giving full information about the deviation from purely periodic sub-pulses.
In Sect. 3 we saw that the components arising in the
2DFS due to the sub-pulse drifting take the form of the Fourier
transform of the combined modulation envelope
,
shifted to
,
.
To obtain the
envelope from the spectrum, we begin by masking out the steady
component
and the mirror-image component. The inverse
transform of this spectrum contains a complex version of the sub-pulse
component of the data which, in the limit of perfect masking
of other
components, is its analytic signal. Its amplitude is only half of
that present in the original data due to the loss of the mirror-image
component, so to circumvent the effect of this in further analysis we
multiply the complex co-efficients of the spectrum by two after
performing the masking. Before performing the inverse transform we
shift the spectrum by the appropriate amount horizontally and
vertically to place the nominal centroid of the sub-pulse response at
DC. The inverse Fourier transform of this shifted spectrum gives the
two-dimensional modulation envelope
,
as used in
Sect. 2.4.
In performing the described procedure, is vital that the mirror image sub-pulse component and the steady component are both carefully removed from the spectrum. In the case where the conjugate mirror image components are not separated by regions of nearly zero power, it is likely that the Fourier transform of the "true'' complex signal had power in opposite quadrants, which as a result of the real sampling is added to that of its conjugated mirror image in the computed spectrum. Under such circumstances, it is not possible to reconstruct the true complex signal of the sub-pulse component, and the derived envelope will differ from that which we had hoped for to a degree dependent on how significant the overlap is. The same caveat applies if the drifting component is not clearly separated from the steady component.
The two dimensional envelope produced using the above procedure can
be decomposed into longitude- and time-dependent parts under the
assumptions made earlier about the form of drifting sub-pulses, as
such:
![]() |
(17) |
With noiseless data and adequate removal of unwanted
spectral components, a convergent result of this algorithm clearly
gives the correct answer. Under the presence of noise, however, care
must be taken in the interpretation of derived amplitude values.
Denoting the rms noise in the real and imaginary parts of the
two-dimensional complex envelope by
(which may be
calculated from the variance of the off-pulse samples in the original
data), the variance of the noise in the real and imaginary parts of
is
![]() |
(18) |
Taking the amplitude envelope incurs some bias due to the presence of
this noise. The derived power values are over-estimated by, on
average,
,
leading to significant bias in parts of
the envelope with low (or zero) signal-to-noise ratio. This bias is
intrinsic to the amplitude-phase representation of noisy data, and so
is also familiar from the estimation of the longitude-dependent
envelope from LRF spectra (Backer 1973) and of modulation indices
(e.g. Taylor & Huguenin 1971).
A second-order effect arises from the fact that
the iterative algorithm described above includes the noise bias in its
unity normalization of the time-dependent envelope, leading to less
correlation and hence an unintended attenuation in the next estimate
of the longitude envelope. The latter effect is easily accounted for
by correction of the normalization to exclude the noise
contribution. Compensating for the former effect requires a more
careful approach.
At the cost of increasing the variance of the noise, the actual
(noisy) values of an envelope may used to determine the bias to
subtract from each element. Alternatively one may take
as an unbiased estimator of the power envelope, however a problem
remains when amplitudes are to be calculated from negative power
values (such as occurs when the signal-to-noise ratio is
low). Nevertheless, this is the approach that is usually employed in
the calculation of modulation indices. A final option that may be of
practical use is to make a model longitude-dependent phase envelope
(since in most cases the measured phase envelope will be
more smoothly varying than the amplitude envelope). Neglecting errors
in the model envelope, the real part of the product
(corresponding in the complex plane to the
projection of the measured envelope on an axis aligned with the model
phase for that longitude) should be free of bias
.
Once the decomposition is performed, four major classes of questions can be asked:
The most popular conceptual framework for understanding drifting sub-pulses is that of a carousel-like system of "sparks'' uniformly spaced about a ring of constant magnetic elevation and altitude from the star surface, with the system as a whole slowly rotating about the magnetic axis (Ruderman 1972). As the pulsar rotates and its beam periodically intersects the observer, a series of pulses, each composed of a series of sub-pulses is seen.
By treating the problem as if all radiation originates on the surface
of a "polar cap'', and is beamed in a direction directly opposite to
that of the center of the star (as shown in Fig. 1),
we may consider the rotational geometry to effect, over the course of
each successive pulse, a sampling of the polar cap emissivity along a
locus defined by the intersection of the sight line and the cap. This
simplification gives a means for good conceptual understanding of the
origin of the observed pulse morphologies and is widely used for this
purpose (e.g. Lyne & Manchester 1988; Mitra & Deshpande 1999; Deshpande & Rankin 2001), but it must be borne in mind
that in reality the radiation is expected to be beamed along a tangent
to the magnetic field line at the emission location. Therefore the
observed emission at a given instant relates to a plasma flux tube
with its foot at some point on an arc on the polar cap connecting the
depicted "line of sight'' and the magnetic pole (labelled
in
Fig. 1c). However, since the field line tangent angle
at a given altitude (where the emission originates) is approximately
proportional to its magnetic colatitude angle, the derived sight-line
locus differs from the true polar cap sampling by a simple scaling of
the magnetic colatitude co-ordinate (see also
Lyne & Manchester 1988)
. In any case, the true
sampling effected in magnetic azimuth is the same as that
derived assuming radial beaming, and if the polar cap emission pattern
is a ring of sparks (or radial "spokes''), it is this co-ordinate
alone that determines the pulse longitude of the observed sub-pulses.
In systems where the sight-line makes a tangential pass of the ring, the emission from the sparks is seen as a sequence of one to a few almost equally-spaced sub-pulses, the phase of which "drifts'' in time with respect to a fiducial point in pulse longitude due to the rotation of the carousel. The number of sparks present, the viewing geometry and the angle between the spin and magnetic axes determine the longitudinal spacing of the sub-pulses (P2), whilst the time spacing (P3) depends on the rotation rate of the carousel and the number of sparks it contains.
There are two main effects predicted under the carousel model that are of importance here.
Firstly, the carousel rotates continuously (with some period
P4 =
NP3, where there are N sparks), meaning that any variations in
spark intensity, shape or spacing about the ring should manifest as a
periodic feature in the drifting sub-pulse signal. This periodicity
would appear as a periodicity in the time-dependent amplitude and
phase envelopes, giving rise in the 2DFS to a pair of "sidebands'' at
frequencies 1/P4 higher and lower than the primary feature at
.
This was the startling finding of
Deshpande & Rankin (2001) who applied the HRFS for the first time, to observations
of PSR B0943+10. The detectability of this effect using 2DFS
techniques is discussed further in Sect. 4.3.3.
The second effect predicted under the carousel model is that
the longitudinal spacing of the sub-pulses depends on pulse longitude.
The sparks are presumed to be spaced uniformly
in magnetic azimuth ()
and their separation in spin longitude
is known, but the mapping from pulse longitude to
magnetic azimuth for the sight line is dependent on the degree of
spin-magnetic pole misalignment (
)
and the viewing geometry
(where
is the angle made between the line of sight and the
positive spin pole). The relation is:
![]() |
Figure 1:
Diagram of sub-pulse emission geometry. In part a) the
angular momentum, magnetic moment, and line-of-sight vectors are shown
with symbols L, ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
Open with DEXTER |
The mapping of pulse longitude to magnetic azimuth is monotonic in the
region in which emission is seen, with a gradient of the opposite sign
to
(where
).
In our analysis, P2 and hence the rate of change of
sub-pulse phase
(
)
are always positive, so
a factor of
will appear in the expression
relating
to magnetic azimuth (
).
The fact that the carousel itself rotates with respect to a fixed
point of magnetic azimuth must also be taken into account. The rotation
period of the carousel is
(see Fig. 1)
where positive
indicates counter-clockwise rotation viewed
from above the active magnetic pole. The direction of drift seen
in the sub-pulses depends also on
:
if their signs are
opposite (as is the case in Fig. 1) the sub-pulses
will drift toward the trailing edge of the profile. We label this
with a negative value for P3, via
.
Since the sub-pulse phase at a given longitude is only sampled once
per rotation period (P1), we will always observe some periodicity
.
Under the carousel model it is possible
that each spark drifts in longitude by more than half of P2 from
one pulse to the next. In this case we may say that the observed
is an "aliased'' version of the true time taken for the
carousel to drift by 1/N turns, P3. The ambiguity in the
"true'' P3 can be expressed in terms of the (signed, nearest
whole) number of sparks (n) that drift by undetected from one pulse
to the next:
.
As shown by Deshpande & Rankin (2001),
this aliasing can potentially be constrained with a measurement of
P4 by observing that
P4/P1 = |N/(n+P1/P3)| and solving for
integer pairs of n and N (in their case
implying n=1 and
N=20)
.
The expression for
is the sum of three terms, the first
a sign-corrected and scaled version of
(Eq. (19)), the second a term which accounts for carousel
rotation over the course of the pulse, and the third an offset
representing the carousel position at
.
![]() |
(21) |
![]() |
Figure 2:
Illustration of the variation in sub-pulse phase (![]() ![]() ![]() ![]() ![]() |
Open with DEXTER |
The techniques used in this work decompose the sub-pulse signal into
the product of two complex envelopes which depend only on time and
longitude respectively. As such, it is important that the longitudinal
spacing of the sub-pulses does not vary with time, and that the
time-spacing of sub-pulses does not vary with longitude. The
longitude-independence of P3 is trivially true in the carousel
model, since the carousel rotates as a whole and the angular speed is
constant around the ring. This is seen in LRF spectra of virtually all
pulsars with drifting sub-pulses, see e.g. Backer (1973). The
time-independence of P2 is technically not satisfied under the
carousel model, since the longitudinal spacing depends on P3 (see
Eq. (20)), which in turn depends on the circulation rate,
which may vary with time. However, for this effect to be significant
Eq. (22) must be dominated by the carousel rotation,
i.e. .
In this case, even minor variations in the rotation
rate could cause the observed (highly aliased)
to vary
wildly in sign and magnitude. To our knowledge, the only pulsar to
clearly show sense reversals in its drift is PSR B0826-34, and the
longitudinal spacing is (within errors) independent of the
instantaneous P3 (Biggs et al. 1985). Other pulsars that
show major changes (but not sign reversals) in the drift rate are
PSR B2016+28 (Taylor et al. 1975) and PSR B0809+74 (after nulls;
Page 1973; Lyne & Ashworth 1983), both of which also show constant longitudinal
spacing. It therefore appears safe to assume that most systems are not
highly aliased and the longitudinal spacing can be assumed constant
with time.
The longitudinal dependence of the observed position angle of linearly
polarized radiation is expected to be intimately related to carousel
rotation sub-pulse behaviour, and comparisons between polarimetric and
drifting-sub pulse analyses have the potential to provide more
information than either alone. Under the "magnetic pole'' or
"rotating vector model'' (RVM; Radhakrishnan & Cooke 1969), the position angle
is that of the sky projection of the magnetic field line on which the
emission region at a given instant lies. Analogous spherical geometry
to that used in the derivation of the magnetic azimuth of the sight line
leads to a similar relation:
In principle,
and
can be determined directly by
fitting the observed position angle swing. These values could then be
used with a measurement of the longitude-dependent sub-pulse phase
envelope and Eq. (20) to determine N and check the
validity of the carousel model (both through its ability to fit the
data and by unambiguous consistency of the fitted value of N with an
integer). In practice there may often be little measurable curvature
in either function, in which case it makes more sense to work with the
(fitted and presumed constant) gradients of the curves.
The first thing to notice is that the position angle swing may have
either positive or negative slope, depending on the sign of .
The maximum rate of position angle swing occurs at the closest point
to the magnetic pole, where
.
If the degree of sub-pulse aliasing is known (or
presumed zero) we can determine the physical direction of
carousel circulation via
.
Secondly, the difference in the absolute value of the slopes can
provide useful information. Since
(for small
), one may take
We note that the work of Deshpande & Rankin (2001) arrived at
, seemingly under the
assumption of
.
Their result of N=20 for
PSR B0943+10 remains valid under the analysis of the preceding
paragraph. Taking their value of 34 for the maximum sub-pulse phase
slope, using the position angle slope of -2.5 from Suleymanova et al. (1988) and
neglecting the (probably significant) uncertainty in this measurement,
we find
.
The marginal consistency may indicate that
is close to -1, implying that
is
close to
and the neutron star rotates clockwise as viewed
by the observer, in agreement with Deshpande & Rankin (2001). If measurements of the position angle and sub-pulse phase slopes
with quoted uncertainties were available, it would be possible to use
the N=20 result
to constrain
and
without the usual recourse to the use
of pulse width-period relations.
Whilst the techniques outlined in Sect. 4.1 are simple and quite direct, we feel that it is important to demonstrate that they perform as expected. To this end, we used the carousel model to produce several simulated time series which we subjected to 2DFS modulation window estimation.
For each sample, we computed the pulse longitude (), and the
corresponding magnetic azimuth (
)
that is sampled by the
sight-line at that longitude using Eq. (19) and a
geometry specified by
and
.
This was
converted to a "sub-pulse phase'' (
)
through multiplication
by
,
and
addition of an offset
(recalling that tis taken as constant in any given pulse), advancing at a rate of one
turn per interval P3(t), which produces the drifting as a function
of time. To investigate the detectability of variations in P3, we
simulated the effect of a pulse null as a period of no drifting
(
)
during the null, followed by an exponential
recovery to the asymptotic rate P3'. The full form of
P3(t) is given by the following equation:
![]() |
(26) |
For each sample the flux contribution from each spark was computed as
a Gaussian in magnetic azimuth (of width
and a peak value
of unity) and the contributions added. The result was scaled by a
factor of
(
)
to simulate an asymmetry
in the azimuthal flux distribution of the spark carousel, to test for
detectability of the carousel circulation time.
Longitudinal amplitude windowing was applied through multiplication by a
Gaussian in
(of width
,
centered
at
with unity peak), and
the resulting intensity was in turn scaled by a factor
![]() |
(27) |
This model leads to a sequence of pulses, each composed of sub-pulses
whose amplitudes vary (in the specific simulations below) from 10% to
of the unity-peaked overall longitude envelope.
The modulation envelopes measured from the simulated results are
compared with a predicted envelope with amplitudes given by half the
peak-to-peak sub-pulse variation produced by the noiseless model at a
given longitude, after normalizing for a unity-mean time-dependent
amplitude envelope. The predicted longitude envelope peaks at an
amplitude of
0.45, compared to
0.50 for the average
profile: that is, a small non-drifting component is present, since the
model carousel has some emission at all points on the ring.
We chose a set of model parameters that give a pulse sequence similar
to those observed from PSR B0809+74 at frequencies around 400 MHz. Specifically, we used a geometry of
,
(Rankin 1993a,b), requiring (from
Eq. (22)) N=16 to reproduce the reported P2/P1 of
40 ms/1.29 s
.
A value of 11.0 P1 was used for P3'.
The main pulse window had a full width at half-maximum (FWHM) of
(from
turns), with each sub-pulse
having a FWHM of
(from
). With the
exception of
and
,
all parameters above derive from
Taylor et al. (1975). We included a null at
tn = 100 P1, of duration
d=8 P1 and recovery time constant
.
No spark-to-spark
intensity variation was included in this model (i.e. c=0) since a
time-varying drift-rate complicates the analysis of this
behaviour
. We produced two data sets, one
with no noise (
)
and one with moderate noise levels
(
), each with 512 pulses in 4096 longitude bins, of
which only the central 512 were used.
![]() |
Figure 3: Simulated PSR B0809+74-like data. A sequence of 100 pulses from the noiseless data are shown, with the interval chosen to show the null. |
Open with DEXTER |
![]() |
Figure 4: Two-dimensional fluctuation power spectrum of simulated PSR B0809+74-like data. To enhance low-level structure, grey-levels have been over-saturated such that pure black corresponds to 10-3 times the peak power. The resolution used here is higher than that of the direct DFT and is achieved by zero-padding the input data. |
Open with DEXTER |
The noise-free data and the resultant 2DF power spectrum are shown in
Figs. 3 and 4 respectively. The form
of each is as expected from the considerations of the preceding
sections. In particular, the spectrum consists of a DC component, and
the first and second harmonics of the drifting response, with all
components significantly extended in the horizontal axis (due to the
longitudinal windowing) and also at a low level in the vertical axis
(due to the nulling-induced time windowing). For clearer
presentation, the spectrum in Fig. 4 was produced
with enhanced resolution by padding the input data with zeroes before
performing the DFT. In practice this would usually not be done, as it
is unnecessary for the determination of a nominal
and
,
and complicates evaluation of the significance of
concentrations of power. We therefore take the position of the
coordinates of the peak bin of the 2DFS computed at the raw
(non-padded) resolution as our nominal periodicity, giving
and
.
This is in agreement
with the parameters used to generate the data, given the bin size of
in
.
After forming the full (symmetrical) complex spectrum and shifting it
by the nominal values for
and
,
the responses
of the average profile and the mirror-image of the fundamental were
removed with a notch filter. To reduce the degree of "ringing'' in
the impulse response of the resultant spectrum, the filters
incorporated a gradual transition of the form
The resultant spectrum was then scaled by a factor of two to account for the use of only one of the two (conjugate-pair) components over which the modulation power is spread, and the inverse Fourier transform was taken. Using the scheme described in Sect. 4.1, we then decomposed the resultant two-dimensional complex modulation envelope into the product of longitude and time envelopes. The results are shown in Figs. 5 and 6.
![]() |
Figure 5:
Complex time-dependent modulation envelope and
residuals for PSR B0809+74-like simulated data. The central panels c)
and d) show the amplitude and phase response expected from the
simulation parameters and the chosen nominal ![]() ![]() |
Open with DEXTER |
![]() |
Figure 6: Complex longitude-dependent modulation envelope and residuals for PSR B0809+74-like simulated data. See caption for Fig. 5 for descriptions of the data shown in each panel. |
Open with DEXTER |
It is evident from the figures that some systematics are present in the residuals of envelopes of both axes. Beginning with the time-dependent envelope of the noiseless data, we see that there is some low-level periodic variation in both the amplitude and phase terms (panels b and e), with a period equal to P3. We believe that this is due to a combination of two effects. The first is the presence of the second harmonic of the sub-pulse modulation, which is shifted in the analysis to the former position of the fundamental. The second is the impulse response of the notch filters used to removed the DC component. This has a similar form to a sinc function (with a somewhat shorter time response due to the choice of transfer function), producing the enhanced oscillations around the null and the boundaries of the data set. The effect of the latter can be reduced by making a smoother transition in the filter, however this necessitates the loss of a greater range of the spectrum, which has the effect of corrupting both the amplitude and phase response to the post-null frequency recovery. In our view the phase evolution is more likely to be of interest than highly accurate amplitude evolution information, so a narrow filter was chosen.
As can be seen in panel (e), there is very little error in the estimated phase, beyond the low-level zero-mean P3 periodicity. The residuals of the noisy data appear, as expected, consistent with the addition of uncorrelated zero-mean noise to the values estimated from the noiseless data. The degree of scatter seen here indicates that estimation noise is likely to dominate over the systematic effects described above for data of all but the highest of signal-to-noise ratios. An exception is the time during the null itself, where the phase is of course not measurable, and the estimated amplitude is affected by the filter response. This error is trivial to correct by examining the pulse-to-pulse mean flux density to identify nulls.
The longitude-dependent modulation window (Fig. 6)
also shows some systematic error in the estimates made from noiseless
data. We believe this is largely due to errors in the estimated
time-dependent phase envelope discussed earlier. This is a
second-order effect, and its magnitude is very small, in the amplitude
window appearing at the one percent level. The systematic trend in
the amplitude window estimated from the noisy data is a simple
artifact of the representation of noisy complex data as amplitude and
phase values, as discussed in Sect. 4.1. Since the variance
of the noise in the synthesized data was 0.25, and given the
doubling of spectral coefficients to compensate for filtering out the
mirror image component, the mean power of the noise in the 2D envelope
is expected to be 1.0 per complex coefficient, giving
.
This leads to a variance of
in the
real and imaginary parts of the longitude-dependent envelope (since
), and hence a
bias of up to
at the edges of the envelope. Techniques
for removing this bias are discussed in Sect. 4.1. The
second-order effect of underestimation by a normalization error
(Sect. 4.1) is expected to amount to only
1%,
making it insignificant and not visible in the plots.
As an alternative that illustrates the applicability of these techniques to other observed phenomena, we produced simulated data from a model with similar parameters to PSR B0943+10.
We assumed a geometry of
,
(Rankin 1993a,b) with N=20 sparks (Deshpande & Rankin 2001). A profile width
of
FWHM (Gould & Lyne 1998) (giving
turns) was
used, with sparks evenly spaced with widths given by
,
as in Sect. 4.3.2. A value of
P1/0.5355 was used for
P3' (Deshpande & Rankin 2001), with a moderately strong dependence on spark
number (c=0.2) applied to give rise to a measurable signal of
carousel circulation. No nulls were included in the simulation, but
the instantaneous P3 values were made to have some variation in
order to model the observed (large but probably finite) Q(
). To this end we produced samples from an
uncorrelated Gaussian noise process and filtered them in the
frequency domain with the function
,
to produce a slowly varying function which was scaled to have a
range of
.
This was added to to the nominal 1/P3' to
give an instantaneous 1/P3(t). As in Sect. 4.3.2, one
noiseless and one noisy (
)
data set was produced. Each
data set contained 512 pulses in 2048 longitude bins, of which only
the central 512 were used.
![]() |
Figure 7: Simulated PSR B0943+10-like data. The first 100 pulses from the noiseless data are shown. The near odd-even modulation is visible, as is the amplitude modulation due to the simulated carousel rotation. |
Open with DEXTER |
![]() |
Figure 8: Two-dimensional fluctuation power spectrum of simulated PSR B0943+10-like data. See caption of Fig. 4. |
Open with DEXTER |
The noiseless data and the resulting 2DF power spectrum are shown in
Figs. 7 and 8. Notice that the
1/P3 feature of the fundamental (at
0.5355 P1) is aliased to a
frequency of
-0.4645 P1. This reflects the fact that if each
sub-pulse is associated with its closest neighbour in pulse longitude
from the following pulse, its drift appears to progress from earlier
to later longitudes. This approach treats all measured
as coming from the range
.
The conclusion from the work of
Deshpande & Rankin (2001) is that each sub-pulse actually drifts 0.5355 times the
sub-pulse spacing towards earlier longitudes with each
successive pulse. Hence, the "nearest sub-pulse'' (
)
convention makes a sub-pulse counting error of n=+1 (see
Sect. 4.2.1). It is important to note that this is
simply a choice of convention. The two cases are indistinguishable
without additional clues as employed by Deshpande & Rankin (2001).
The 2DFS was filtered with a function of the form shown in
Eq. (28), with
w=0.03/P1, and shifted to place the observed
peak of the fundamental (
,
)
at
DC. This was scaled by a factor of two and inverse transformed and the
result decomposed into time- and longitude-dependent modulation
envelopes as in Sects. 4.1 and 4.3.2. The
results are shown in Figs. 9 and 10.
![]() |
Figure 9: Complex time-dependent modulation envelope and residuals for PSR B0943+10-like simulated data. See caption for Fig. 5. |
Open with DEXTER |
![]() |
Figure 10: Complex longitude-dependent modulation envelope and residuals for PSR B0943+10-like simulated data. See caption for Fig. 5. |
Open with DEXTER |
As with the 0809+74-like data, some systematics are present in the
results from the noiseless data. In this case they are less severe due
to the good separation between the fundamental of the drifting
component and the DC component. In any real (noisy) data the errors
are likely to be dominated by noise. In this model, due to the
broader simulated main pulse, the phase rate (i.e. )
variation
across the profile is noticeable even in the noisy data.
As expected,
the bias in the wings of the longitude-dependent amplitude envelope is
of the same magnitude as that in the 0809+74-like data.
The amplitude variations due to the simulated carousel circulation are
clearly visible in the modulation envelope estimated from the noisy
data, as are the phase variations due to the simulated frequency
noise. Since the phase variations are quite small, the amplitude
modulation appears highly periodic and can be detected in the power
spectrum of the inferred amplitude envelope
(Fig. 11). The modulation also appears as a
pair of "sidebands'' around the steady and drifting components in the
2DFS (Fig. 8). After previous authors
(e.g. Deshpande & Rankin 2001), we have "stacked'' the columns of the power
spectrum corresponding to
1/P2 < 64 to produce an "average''
fluctuation power spectrum.
Figure 12 shows the stacked power spectrum
around the
-0.4645 P1-1 drifting component. The
sidebands associated with the carousel rotation are clearly present,
offset by
from their
respective parent features.
![]() |
Figure 11: Power spectrum of the time-dependent amplitude envelope estimated from the noisy data. Values are normalised by the mean of noise-floor bins. |
Open with DEXTER |
![]() |
Figure 12:
Portion of "stacked'' 2DFS around the drifting component.
Values are normalised by the mean of noise-floor bins. The central feature
peaks at a value of ![]() |
Open with DEXTER |
The significance of detections of the carousel circulation measured
via the power spectrum of the inferred amplitude envelope versus the
stacked 2DF power spectrum must be evaluated carefully. Since the
former is estimated from a single Fourier transform, the baseline
noise power (suitably scaled) is
distributed with
two degrees of freedom. In contrast, the stacked spectrum is formed
by adding 16 columns of the two-dimensional spectrum, so the (suitably scaled) spectral values follow a
distribution with
32 degrees of freedom. For a detection of 97.8% confidence
(corresponding to the 2-sigma point if the noise was Gaussian), the
threshold normalised spectral power values are
3.8 for the
power spectrum of the inferred amplitude envelope, versus
1.6for the normalized stacked spectrum. The corresponding "4-sigma''
points are
10.4 and
2.3. Neglecting effects due to phase
noise, each responds to a modulation of amplitude c with a spectral
power value proportional to c2, plus the "noise floor'' value of
(on average) unity. From the values of peaks observed in these
simulated data we can infer a proportionality constant of
20between the (noise-subtracted) responses of the two types of
spectra. Hence, values of c giving rise to "2-'' or "4-sigma''
detections in the stacked spectrum (with power levels of 1.6 and 2.3)
would produce far more significant detections (with power levels of
13 and
27) in the power spectrum of the inferred
amplitude envelope. This enhancement in sensitivity (amounting to a
factor of
2 in the minimum detectable amplitude of carousel
circulation modulation) can be understood as a result of exploiting
the inherent phase relations of the complex spectrum, rather than
(incoherently) summing spectral power values.
As noted earlier, the model assumed here also produces sidebands around the DC component of the spectrum. Such sidebands were not observed in the studies of Deshpande & Rankin (2001), a fact that needs to be understood if carousel circulation is to explain the other sidebands. One possible explanation is that the total (integrated) luminosity of each spark of the carousel is approximately equal (leading to no sidebands around DC) but that the width of sparks differs as a function of position in the carousel.
We have shown that the two-dimensional Fourier spectrum of pulsar longitude-time data is of value in the analysis of drifting sub-pulses. We consider the drifting component of the signal as the convolution of a complex "modulation envelope'' with a pure two-dimensional periodicity. The modulation envelope describes the variations in average sub-pulse amplitude and phase as a function of longitude and time, and can be decomposed into the product of two one-dimensional envelopes which are functions of time and longitude respectively. This makes the technique well-suited to studying a variety of phenomena associated with drifting sub-pulses, including longitudinal variations in sub-pulse spacing (induced due to viewing geometry or other factors), variations in the drift rate as a function of time (due to the recovery from a null, random phase noise, etc.), comparison of the longitudinal amplitude dependence of the drifting and steady components of the pulsar emission, and variations in the average sub-pulse amplitude as a function of time (due to nulling, carousel rotation, etc.).
Acknowledgements
We thank R. Ramachandran and M. van der Klis for helpful comments on the text, A. Deshpande for useful discussion, and the referee (J. Middleditch) for drawing our attention to the sub-pulse modulation of "PSR1987A''. RTE is supported by a NOVA fellowship. BWS is supported by NWO Spinoza grant 08-0 to E. P. J. van den Heuvel.
As noted earlier, the 2DFS and HRFS are in fact equivalent. In this section we briefly show how this fact arises.
Consider the response to a single sinusoid,
.
This will appear as a delta function in the HRFS at
,
(both in units of cycles per
P1-interval), where
and
denote the
fractional and integer parts of their arguments. For the 2DFS, when
the data are stacked the signal will appear in each pulse period (of
length P1) as a sinusoid with frequency
(cycles per
radian of pulse longitude) and phase
where t/P1 is the integer pulse number. Transforming first
across rows results in a delta function at
with a phase angle of
(and a corresponding component at
which can be
ignored for now as it is simply the source of the symmetry of the
2DFS). Considering the column corresponding to
,
the response is a pure complex exponential with a frequency of
(cycles per time unit). Therefore the result of
the two-dimensional DFT is a delta function at
,
.
A 1:1 mapping of
,
thus applies in this case, when the fact that the
bin size in the discrete 2DFS is
is considered. Since
both transformations are linear and 1:1 invertible, we therefore
conclude that they are identical in all cases, given appropriate
mapping between parameters.
The derivation of Eq. (19) follows a straightforward
coordinate frame transformation, from the spin frame (in which the
sight line parameters
and
represent longitude and
co-latitude), to the magnetic frame. In each frame longitude is
defined such that the positive spin pole of the complementary frame is
at longitude zero. We definine basis vectors
and
for the
spin and magnetic frames respectively, aligning
(and
)
with the positive poles and placing
and
at longitude zero in their respective frames. In what
follows all vectors (including
and
)
are taken as
being of unit length. The expression for
in the spin frame
is:
![]() |
(B.1) |
![]() |
= | ![]() |
|
= | ![]() |
(B.2) |
The derivation of the polarization position angle follows a similar
path. It is based in a third frame which is fixed to local rest frame
the observer, has its origin at the star center, positive pole
pointing to the viewer (i.e.
)
and places the
positive spin pole at longitude zero. All meridians in this frame
project as lines on the sky with position angles offset from the
position angle of the spin axis by their longitude (
).
Under the RVM, the position angle of
linear polarization is given by the projection of the local magnetic
field line on the sky, which is given by the projection of the great
circle co-planar with the
and
.
Hence the
(counter-clockwise) polarization position angle (minus the position
angle of the projected spin axis) is given by the longitude of the
magnetic pole in this frame. The transformation from the spin frame
is a rotation of
in longitude, followed by a rotation about
by
(versus
in the above) and a rotation of
radians about
as above, whilst the input vector has
spherical coordinates
(versus
in the
above). From this it follows that Eqs. (19) and (23) differ only through
switching of
and
and a sign reversal (in
or equivalently, in the final longitude value).