© ESO, 2012

1. Introduction

As radio emission passes through the interstellar medium (ISM) it interacts with electrons, which cause it to be dispersed. This is observable through pulsed emission, where a low-frequency pulse is delayed with respect to the same pulse at higher frequencies. For a pulse at frequency ν (frequency will be given in MHz here and throughout the paper) travelling a distance, D, through an unmagnetised ionised gas, the dispersive delay with respect to infinite frequency, Δt_DM, is given by: (1)where v_g is the group velocity: (2)ν_p is the plasma frequency: c is the speed of light, m_e and e are the charge and mass of an electron respectively, and n_e is the electron density in cm^-3. Normally, this time delay can be approximated using the first term in the Taylor expansion of 1/v_g, giving the cold dispersion law¹: (5)where DM is the dispersion measure, the integrated column density of free electrons in pc cm^-3. This relation was used as early as 1968 to accurately predict the dispersive time delay to within 1 part in 3000 between 40 MHz and 430 MHz (Tanenbaum et al. 1968). However, as ever higher timing precision is required for projects like using pulsars for gravitational wave detection (Jenet et al. 2005) it is important that any second-order effects of the ISM, such as refractive delays, DM variations and delays associated with pulse broadening from scattering (Foster & Cordes 1990; Cordes & Shannon 2010), are also studied and understood fully (You et al. 2008; Hemberger & Stinebring 2008). Most of these proposed effects have strong frequency dependencies, with scaling indices between ν^-3 and ν^-4, and are therefore most prominent at low frequencies.

As well as second-order ISM effects, there are other proposed frequency-dependent effects, such as propagation from within the pulsar magnetosphere (Michel 1991), super-dispersion (Shitov & Malofeev 1985; Kuz’min 1986; Shitov et al. 1988; Kuz’min et al. 2008) and aberration and retardation (Cordes 1978). The ionosphere also contributes slightly to the DM but typically the total number of electrons along a path through the ionosphere is less than 100 TEC units (e.g. Liu et al. 2009), which corresponds to a DM of just 3.24    ×    10^-5 pc cm^-3. In practice, this contribution is small and follows a ν^-2 law, making it indistinguishable from the dispersive delay of the ISM.

As the sizes of these effects have not yet been determined, their frequency dependence could introduce systematic errors when data from different frequencies are combined. The size of these effects could also potentially vary with time. For example, if the distribution of free electrons along the line-of-sight changes, the magnitude of some of these effects will vary, introducing further errors into pulsar timing data.

Another potential source of error arises from pulse profile evolution. Most (if not all) pulsars show variation in the shape of their pulse profile as a function of frequency (Craft 1970). This change can be anything from a slight broadening of the pulse, or difference in the relative positions of the components, to components appearing or disappearing completely.

Even though wide-band simultaneous observations have been possible to perform for a long time (see for example Phillips & Wolszczan 1992; Kuz’min et al. 1998; Hankins & Rankin 2010) the relatively narrow frequency bands compared to the separation interval between frequencies have meant identifying exactly how the components evolve has been difficult, particularly at the lowest frequencies. One of the reasons for this is that profile evolution, which only manifests itself clearly over very wide frequency ranges, is difficult to disentangle from dispersion (dispersion is normally removed by fitting the data with a ν^-2 power law which mimics Eq. (5)). For example, in Kuz’min et al. (1998) the authors suggest that PSR B0809+74 appears to have an extra, non-dispersive delay of  ~ 30 ms at 10 GHz, though this could alternatively be explained by pulse profile evolution.

To calculate pulse times of arrival (TOAs), the data are cross-correlated with a “template”, a smoothed model of the pulse profile. The peak of this cross-correlation spectrum is used to determine when the pulse arrives relative to the template, and the phase of the Fourier transform of the cross-correlation is then used to measure the phase offset. If the pulse profile evolves significantly with frequency the shape of the data and the template can be slightly different. These subtle changes in shape can introduce frequency-dependent errors into TOAs if they are not accounted for properly. Ahuja et al. (2007) simulated the effects of pulse profile evolution on measurements of DM. They found that using a template which is slightly different from the shape of the data can cause a gradient in the phase of the cross-correlation, which causes an offset in the TOA.

Some of the simulated errors, in the work of Ahuja et al., due to this effect in normal (slow) pulsars were as large as milliseconds, although the magnitude of the effect on real data has not been studied in detail until the work presented here. The size of this effect should scale linearly with pulse period, and so it is expected to be much smaller for millisecond pulsars, which is why sub-microsecond timing precision is already common (for example, see Verbiest et al. 2009). But, with the introduction of the first ultra-broadband receivers to be used for pulsar timing (e.g. DuPlain et al. 2008; Jones et al. 2010), which will be capable of observing with large fractional bandwidths even at frequencies of a few gigahertz, pulse shape variation across the band will become far more pronounced, and this effect will become more apparent.

We took observations of four bright pulsars (pulsars B0329+54, B0809+74, B1133+16 and B1919+21) simultaneously at multiple frequencies with the LOw Frequency ARray (LOFAR, see van Haarlem et al., in prep.), the 76-m Lovell Telescope and the Effelsberg 100-m Telescope, spanning a frequency range between 40 and 8350 MHz to try to constrain the properties of the ISM and test the cold plasma dispersion law. Details of these observations are given in Sect. 2. Our observations cover a very large frequency range simultaneously and each recorded band has a significant fractional bandwidth, also making them useful for studying pulse profile evolution. LOFAR provides the largest fractional bandwidth which has ever been possible at the lowest radio frequencies observable from Earth.

The observations from both LOFAR observing bands were also taken truly simultaneously, both starting and ending at exactly the same time², and passing through the same hardware so no timing model was required to align the profiles precisely (as opposed to using overlapping observations with a timing model). This part of our frequency range provides the best test of steeply scaling frequency-dependent delays and this allows us to place strong constraints on the magnitude of the second-order ISM effects. Any deviation from the ν^-2 law larger than a few milliseconds would be clearly visible over this large bandwidth and at these low frequencies.

Similar observations have been carried out previously by Tanenbaum et al. (1968), who were able to set a surprisingly good limit that the error in the dispersion equation, Δt < 3    ×    10^-4 s between 40 MHz and 430 MHz just six months after the discovery of pulsars in 1967. Since then similar work has been done by Shitov & Malofeev (1985), who found evidence for excess dispersion at low frequencies. This was also supported by Hankins et al. (1991), who noticed a difference in the DM values determined from aligning single pulses and those determined from aligning the average pulse profile. These delays at low frequencies have been termed “super-dispersion” and attributed to magnetic sweepback in the pulsar magnetosphere (Shitov 1983). Subsequent work by Phillips & Wolszczan (1992) and Ahuja et al. (2005) found contradictory results, with no evidence for extra dispersion at low frequencies, but some pulsars (B0525+21, B1642-03 and B1237+25) showed an increased dispersion measure at high frequencies. These high frequency delays were explained by propagation delays in the pulsar magnetosphere.

In this paper, after detailing the observations (Sect. 2), our analysis (Sect. 3), and the simulations used to determine upper limits on the magnitude of the delay present in our data (Sect. 4); we determine the effect of the ISM on pulsar timing and extract information about the composition of the ISM (Sect. 5) and the pulsar magnetosphere (Sect. 6). We also use these data to determine how much of an impact profile evolution has on pulsar timing and use the pulse profiles in each of the frequency bands to construct a model of how the pulse profile evolves with frequency (Sect. 7). Finally, we summarise our main findings in Sect. 8.

2. Observations

On 11 December 2009 we observed four pulsars using LOFAR, the 76-m Lovell Telescope and the Effelsberg 100-m Telescope simultaneously. LOFAR is an international interferometric telescope, comprised of many thousands of dipole antennas grouped into “stations” and operating in the lowest four octaves of the “radio window” visible from the Earth’s surface. The LOFAR stations are arranged in a sparse array, spread across Europe, with a dense core region located in the Netherlands. LOFAR operates in two frequency bands which straddle the FM radio band. The Low Band Antennas (LBAs) can record 48 MHz of bandwidth between 10 and 90 MHz and the High Band Antennas (HBAs) can record 48 MHz of bandwidth between 110 and 240 MHz. For a full description of how LOFAR is used for pulsar observations see Stappers et al. (2011) and for a general LOFAR description see van Haarlem et al. (in prep.).

At the time of the observations, LOFAR was still under construction so only single stations were used to record data and the LBAs were only able to record 36 MHz of bandwidth. A core station (CS302, which consists of 48 HBA tiles) was used to observe in the LOFAR high band between 138.9709 and 187.4084 MHz and an international station (DE601, based at the Effelsberg Radio Observatory and consisting of 96 LBAs) was used to observe between 42.0959 and 78.4119 MHz in the low band. It is worth noting that the LOFAR stations which were used were not yet internally calibrated and hence the sensitivity of the observations presented here is significantly lower compared to what is now possible.

Data were taken simultaneously with the 76-m Lovell Telescope, using the digital filterbank backend in search mode, at a central frequency of 1524 MHz with 512 MHz bandwidth (see Hobbs et al. 2004, for details). Simultaneous data were also taken with the Effelsberg 100-m Telescope, using the Effelsberg-Berkeley Pulsar Processor to coherently dedisperse the data online, at a central frequency of 8350 MHz with 1000 MHz bandwidth (see Backer et al. 1997, for details). The observational parameters are summarised in Tables 1 and 2.

In this paper, we concentrate mainly on the total intensity profiles of pulsars, and only Stokes I data were recorded with each instrument (LOFAR, the Lovell telescope and the Effelsberg telescope). However, if there is a significant gain difference between the two hands of (either linear or circular) polarisation, or leakage, and the pulsar is strongly polarised then there may be distortions in the profile. Instrumental polarisation from the Lovell and Effelsberg telescopes is known to be small (see e.g. Montenegro-Montes et al. 2008; Gould & Lyne 1998), but at the time of these observations, the LOFAR polarisation data was uncalibrated. We note however that in the case of LOFAR the leakage terms are at worst a few percent, and this is further reduced due to the fact that there are many thousands of elements, which are physically identical, so any leakages would average out over the array. Also, as the orientation of the LOFAR dipoles is at 45 degrees to the north-south orientation and our sources were all observed close to transit, both sets of dipoles would have received approximately the same amount of radiation from the source. Lastly, we note that the polarised fraction is relatively low, and that the Faraday rotation in the ISM is large at these frequencies. In general, with the bandwidths used to make pulse profiles (typically  ~ 12 MHz, see Sect. 3), the plane of polarisation is rotated through at least 180° in each profile, thereby further reducing the effect of polarisation calibration terms. We are therefore confident that the results presented here are not affected significantly by polarisation calibration.

3. Analysis

The data were converted into PSRFITS format and processed using the PSRCHIVE software suite (Hotan et al. 2004). Regular timing observations from the Lovell telescope at Jodrell Bank Observatory (Hobbs et al. 2004) were used to derive accurate values for the spin period (P) and spindown (Ṗ) for the day of our wide-band observing campaign. Initial DM estimates were taken from Hobbs et al. (2004, see Table 2. The data were dedispersed, folded and terrestrial interference signals were removed by hand using the RFI-removal software, pazi³.

Templates were made for the observations at each telescope by completely collapsing the data in both time and frequency to make an average pulse profile. The profiles were then fitted with von Mises functions⁴ to create analytic templates using paas³.

A template was created for each observing band and these templates were subsequently aligned by eye with pas³ using either the brightest peak (for PSR B0329+54) or the midpoint between the two brightest components (for pulsars B0809+74, B1133+16 and B1919+21) as the fiducial point as described in Craft (1970). The templates were cross-correlated with the data to get times of arrival (TOAs) using pat³ and barycentred using tempo2⁵. These TOAs were subtracted from a model of the pulsar’s rotation in order to produce “timing residuals”, which are plotted in the left hand panel of Fig. 1.

Currently, it is not possible to perform absolute timing with LOFAR. The clock corrections between the stations are known to better than a few nanoseconds, but the offset between the LOFAR clocks and Coordinated Universal Time (UTC) is not logged. This introduces time delays which are on the order of a few milliseconds between LOFAR data and the data from the other telescopes in the observations. The delays presently require an arbitrary phase offset (“jump”) to be removed between data sets from different telescopes. In the case of the 8.35 GHz Effelsberg observation, where no frequency resolution was stored (see Table 1), the determination of this phase offset uniquely defines the residual with respect to the other observations, implying that no further timing information can be derived from the observation. We have therefore omitted the Effelsberg observations from the timing analysis (though they were still used in the analysis of profile evolution, see Sect. 7).

One can see in Fig. 1 (left) that for each of the pulsars there is clear structure in the timing residuals as a function of frequency. Often, this structure shows different slopes for each of the bands, indicating that the deviations from a good fit are not simply caused by an incorrect dispersion measure. Initially it seemed likely this might be caused by deviations from the simple form of the cold plasma dispersion relation.

Consequently, we tried fitting the data using power laws with exponents between  − 5 and  + 5, but none of the fits significantly improved upon the chi-squared obtained from only fitting for the dispersive delay (for example, in the case of PSR B1919+21, the best power-law fit only decreased the reduced chi-squared from 1.47 to 1.44, which, given the number of degrees of freedom in the model, still corresponds to a  ~ 99% chance that there is unmodelled structure in the residuals). In some cases, the residuals also showed structure which could not be explained by a single function (for example, a positive gradient in the LBA data, a negative gradient in the HBA data and a positive gradient in the data from Jodrell Bank). No simple power law can account for the different slopes in our timing residuals, so we conclude that, the structure in the residuals is not caused by the ISM or aberration and retardation in the pulsar magnetosphere.

The real cause of the systematic errors is the pulse profile changing as a function of frequency, as noted in Ahuja et al. (2007). As the profile changes, the template which is used in the cross-correlation becomes more and more inaccurate and there is a systematic drift in the residuals towards earlier or later TOAs. Figure 2 shows how these errors arise. The top panel shows the template used for the HBA observation of PSR B1919+21 along with the observed pulse profile at 145 MHz. The bottom panel shows the cross-spectrum phase (see Ahuja et al. 2007, for details) of the two profiles and a straight line fitted to the data. There is a gradient in the fitted line which means that the peak of the amplitude of the cross-correlation (used to produce the TOAs) is shifted slightly, causing the apparent TOA to be delayed or advanced compared to its true value.

This effect has only been brought to attention recently by Ahuja et al. (2007), because pulsar timing is normally done at high frequencies (where the pulse shape changes are less obvious), and with modest bandwidths. With these narrow bandwidths, it is easy to incorrectly interpret this structure as an incorrect DM. With the wide fractional bandwidth and high frequency resolution of our simultaneous data, however, the effect cannot be fitted with a simple ν^-2 law, and so the structure is more easily identified.

To remove the effects of profile evolution and find the true DM value and residuals, we first modelled how the pulse profiles evolve as a function of frequency. The data were divided into narrower-frequency subbands. Each subband was chosen such that the profile evolution within it was not significant, but there was still enough bandwidth so that the pulsar was detected clearly. The size of these bands depended on the signal-to-noise ratio of the observation, but typically we used 4 subbands per telescope, corresponding to  ~12 MHz in the LOFAR data and  ~100 MHz in the Lovell data. Each subband was then collapsed in frequency and folded in time to create a pulse profile, which, unlike the previous analysis, we fitted with Gaussians⁶.

The width, height and positions of the Gaussians were free to vary, and the best fit values for each of them were recorded and used to create a continuous model of how the pulse profile evolves as a function of frequency. The fiducial point was chosen as either the centre of the main pulse or the midpoint of two brightest components depending on the morphology of the pulse profile (see Craft 1970).

This process is shown for PSR B1919+21 in Fig. 3. We used power laws to model the frequency-evolution of each of the fitted parameters, as they are the simplest functions which fit the data well. However, as can be seen in Fig. 3 the fits were not perfect. The separation of the peaks follows a power law well (although the exponent is positive, contrary to what is expected by radius-to-frequency mapping, see Sect. 7.2.4), but the widths of the components are not well described by power laws. This is due to the overlap between components. In the overlapping region it is unclear how much power belongs to each component, and the solution which is achieved from the fit is not unique. The ratio of the peaks also could not be fitted with a single power law, which was the case for all of the pulsars in our sample. Instead we used a power law for each frequency band. This is because both of the components have different spectral profiles (which may contain one or more breaks) and as we don’t have absolute flux values, it is unclear whether one component is getting brighter or the other is getting fainter. Scintillation across the different frequency bands could also contribute slightly to the discontinuities in the component amplitudes, although the effects of scintillation will be small in the average pulse profiles because of the long integration times used in these observations.

Our model of the pulsar was then used to create a template for each frequency channel, which was cross-correlated with the data to get TOAs using pat in the same way as the single templates were.

The right hand side of Fig. 1 shows the residuals obtained from using the templates derived from this frequency-dependent model. Using these model-based templates for timing improves the residuals and reduces the systematic errors which were seen in the initial residuals. With the frequency-dependent templates the residuals show far less systematic trends and agree to within the error bars on a single DM in most of our observations (see Fig. 4 and Table 3). PSR B1133+16 seems to show some structure in its residuals. Its DM at  ~1400 MHz is more than 3 sigma from its DM at low frequencies and this cannot be removed by using a frequency-dependent template. This is possibly due to the finite sampling of the data and template which means that at the highest frequencies the pulse profile of PSR B1133+16 (whose components are the narrowest of all the pulsars in our sample) is barely resolved. Subsequent observations taken at Jodrell Bank, with higher time resolution, show no such structure. None of the other pulsars in our sample were affected by this problem.

We will return to the issue of the component evolution in pulsars B0329+54, B0809+74, B1133+16 and B1919+21 and its impact on pulsar timing (Sect. 7), but first we will consider the influence of the ISM.

4. Simulations of non-dispersive delays

Before we started looking for extra structure, the data had already been fitted for DM and a phase offset between LOFAR and the Lovell telescope. The size of the error bars from the cross-correlation are also very different in each of our observing bands. For these reasons, it was necessary to perform simulations to determine how sensitive the data truly are to these frequency-dependent effects.

To simulate the effects of the ISM, we added a delay to the folded profile in each subband of our data according to a ν^-4 scaling law (most second-order ISM delays scale roughly as ν^-4, see Sect. 5). This data was then cross-correlated with a frequency-dependent template, fitted for dispersion measure and a jump between the LOFAR data and the L-band data, and displayed in the same way as the real data to produce residuals which had an artificial delay added to them.

We used the least squares method to fit the residuals with a ν^-4 power law. The chi-squared value of the fit and the number of degrees of freedom were computed, and used to compare the likelihood that the structure in the residuals was caused by a ν^-4 delay with the likelihood that it was caused by chance⁷. We reduced the magnitude of the delay which was introduced to the data until the residuals had an equal probability of being caused by chance as being caused by an ISM-like power law. This delay was set as the upper limit on the magnitude of this effect in our timing residuals. This process is demonstrated in the left panel of Fig. 5. From these simulations, we determined that our sensitivity to steep, negative, frequency-dependent power laws is dominated by the error bars associated with the LBA observations.

In the same way, we also simulated the effects of an aberration/retardation-like delay. Cordes (1978) showed that the emission height, r(ν) ∝ ν ^− 2x, where x is a constant which depends on the radius-to-frequency mapping of the particular pulsar. The author then used data to show that typically x falls in the range 0.21 − 0.55.

The delay in the pulse arrival time is proportional to the height at which it is emitted (see Eq. (19)). From this, we can deduce that Δt should scale with ν^0.4 − 1. For our simulations, we elected to use the average value of x found by Cordes, so chose a power law which scales with ν^0.6.

Again, we compared the likelihood of a fitted function and no function, and reduced the magnitude of the effect until it was undetectable in our timing residuals (see the right panel of Fig. 5). Remarkably, the fitted jump has very little effect on our ability to detect aberration and retardation effects, and because the error bars in our high frequency data are much smaller than those at low frequencies we are, in fact, much more sensitive to high frequency delays. If the power law is steeper than ν^0.6 our sensitivity to these effects increases. Table 4 shows the determined upper limits on high and low frequency delays in our data.

5. ISM effects

5.1. Impact on pulsar timing

There are a number of second-order effects caused by the ISM which could induce additional delays in our data. These effects scale strongly with wavelength and are at their strongest at low frequencies. As shown by our simulations, the lack of any remaining structure in the residuals indicative of unmodelled effects shows there is no significant indication of super-dispersion, refraction, anomalous dispersion or frequency-dependent DM variations. We can, however, still use these observations to place important limits on the magnitude of some of these effects, at least for these four sources.

Taking the maximum unmodelled (i.e. non-dispersive) ISM delay in the LBAs and extrapolating it to higher frequencies gives an indication of the impact of the ISM on pulsar timing projects like Pulsar Timing Arrays (PTAs, see Romani 1989; Foster & Backer 1990). Figure 6 shows the magnitude of ISM effects as a function of frequency for two pulsars with the largest (B0809+74) and smallest (B1919+21) deviations from white noise in our sample. The largest possible delay in our data (see Table 4), scaled with the gentlest gradient (i.e. the largest possible deviation) still corresponds to only  ~50 ns at 1400 MHz. Although this figure will change significantly for every line-of-sight, depending on the parameters of the ISM (and potentially when the observation was taken) this upper limit is only roughly half of the  ~100 ns precision required for the first generation of PTAs (Jenet et al. 2005).

Cordes & Shannon (2010, and references therein) describe a number of frequency-dependent delays, which are caused by the distribution and number of free electrons along the line-of-sight. Combining these relations with our timing residuals allows us to constrain some properties of the ISM (see Table 5).

5.2. Frequency-dependent DM variations

Frequency-dependent DM variations are caused by scattering by the ISM. The diameter of the scattering disk (the region in the sky where scattered radiation is received from) increases as  ~ν^-2.2. The DM, which is an average over the scattering disk weighted by the flux received from a given direction, therefore also varies as a function of frequency (Fig. 7).

If we assume that all of the scattering material is concentrated in a thin screen half-way along the line-of-sight, the variation in DM is given by (Cordes & Shannon 2010): (6)where D′ is the distance between the pulsar and the scattering material (in kpc), ν is the observing frequency in MHz and SM is the Scattering Measure along the line-of-sight in units of kpc m^−20/3. This corresponds to an RMS variation in time delay of: (7)This allows us to constrain the SM along the line-of-sight (see Col. 4 of Table 5). The SM must be fairly large for this effect to be detectable and all of the pulsars which we observed are at low DM (<30 pc cm^-3) and show very little scattering (with the exception of PSR B0329+54, no scattering is visible in any of the pulse profiles). Our upper limits are about 2–3 orders of magnitude above those from the NE2001 model of Cordes & Lazio (2002).

To detect this effect in PSR B0329+54 (assuming the SM from NE2001) would require a timing residual RMS of 10 μs, so it is not surprising that this effect was not detectable in our data. However, it is possible that this effect could be detected in the future with LOFAR using pulsars with a higher SM (for example PSR J2044+4614 or PSR B2036+53, both of which are closer to the Galactic plane). The LOFAR data presented here is far from full sensitivity, and, when the data are fully calibrated, and all 24 of the stations in the LOFAR core can be combined coherently, the telescope’s sensitivity will be improved by at least an order of magnitude.

5.3. Refraction

Refraction can also introduce a delay into pulsar timing data at low frequencies. If there is a gradient in the interstellar electron density perpendicular to the line-of-sight, rays of light are refracted, and bent through some angle θ_r. For a thin screen this is given by Clegg et al. (1998): (8)where r_e is the classical electron radius and N_e is the column density of electrons along the line-of-sight. At lower frequencies, light is bent through a larger angle. The larger the angle, the further the light must travel and the longer it takes to arrive at the observer (see Fig. 8).

Cordes & Shannon (2010) showed the delay between a straight path and a refracted path is given by: (9)where D_eff is a characteristic distance, given by: (10)D is the distance between the observer and the pulsar and D_s is the distance between the observer and the scattering screen. Assuming the line-of-sight is dominated by one screen (which is a reasonable assumption given the evidence in Stinebring 2006; Brisken et al. 2010, and that our sample is quite nearby), the thin screen model is valid. By substituting θ_r from Eq. (8) into (9) and rearranging for ν we can derive an expression for the time delay between refracted and unrefracted rays: (11)where is the average gradient in the DM perpendicular to the line-of-sight in pc cm^-3 AU^-1, ν is in MHz and D_eff is in kpc . This number can be constrained from the lack of structure in our timing residuals and can be used as an indicator of how much the DM of a pulsar is likely to change.

It follows from the derivation above, that the value obtained holds over distances dx of the size of the scattering cone. At these frequencies the scattering cone for a typical (nearby) pulsar is on the order of a few AU (for B0329+54 the scattering cone is 10 AU at half the distance to the pulsar), so the most natural units to use for this quantity are pc cm^-3 AU^-1. The unit pc cm^-3 AU^-1 also corresponds to how much the DM of a pulsar at a distance D travelling 250 km s^-1 will change in approximately one week. It should be noted however, that the number is an average over the entire scattering disk, and so is not sensitive to small scale anisotropies. Upper limits are shown for for the pulsars in our sample in Col. 5 of Table 5.

The main source of error in this number arises from D_eff as the distance to the scattering screen is unknown. The function is very sensitive to nearby scattering screens (D_s < 0.25D) but is not very sensitive to distant scattering screens (D_s > 0.75D). It is, however, approximately constant for 0.25 < D_s < 0.75, varying by only a factor of  ~3 in this range. In our calculations, we assume that the scattering screen is roughly half way along the line-of-sight and D_eff ≈ D.

5.4. Anomalous dispersion

By modifying Eq. (1) to include gyro-motion in a magnetic field and electron collisions, and including the quartic term in the Taylor expansion of 1/v_g, we can write a more general version of the plasma dispersion law, which can be expressed as the sum of three functions (Eqs. (13), (15) and (16), see Phillips & Wolszczan 1992): (12)The γ₁ term: (13)is the normal dispersive delay term multiplied by an additional factor which depends on the ratio of the thermal velocity of the electrons to the speed of light, β_therm: (14)Where k is the Boltzmann constant, T is the temperature of the plasma and m is the mass of an electron. This means that if the temperature is non-zero it adds a small contribution to the dispersive delay. In practice this is a very small number and is absorbed into the ν^-2 law making the DM slightly higher than it would be normally, but not altering the power law of the delay. Even an unrealistic change in temperature as large as 100 K only modifies the apparent DM by 0.01%, and would be indistinguishable from the normal dispersive delay.

The γ₂ term is because of the plasma being weakly magnetized. For a circularly polarised wave: (15)where RM  is the Rotation Measure in rad m^-2. This only affects circularly polarised sources and is generally negligible. The RM along a given line-of-sight is closely related to the SM and the DM, so when a pulsar’s RM is large enough to make this effect detectable, scatter broadening makes the pulsar undetectable. An RM of 280 would be needed to produce a detectable delay of 1 ms at 20 MHz. The scattering time for a pulsar with this RM would be  ~4 s, which, in most cases, would render the pulsar undetectable (assuming RM ~ 5    ×    DM, the average from the ATNF pulsar catalogue Manchester et al. 2005; and the empirical scattering law from Bhat et al. 2004).

The γ₃ term depends on the Emission Measure, EM  along the line-of-sight: (16)where EM is given in pc cm^-6. For a uniform distribution of plasma between the pulsar and observer the contribution from the ν^-4 term is small even at low frequencies and is probably not detectable. However if the electrons are arranged in a clumpier distribution along the line-of-sight (as discussed in Kuz’min et al. 2008) this term becomes larger.

Whilst both the β_therm and γ₂ terms are too small to detect, the fact that no delay is observed in our data can be used to constrain the EM along the line-of-sight through the γ₃ term. The upper limits are given in Col. 6 of Table 5. Although they are currently around five orders of magnitude higher than the values predicted by NE2001 (Cordes & Lazio 2002), using this method on observations at lower frequencies and with higher sensitivity, could provide a new way to measure the EM along a given line-of-sight.

6. Magnetospheric effects

In addition to the delays caused by the ISM, there are potential sources of delay from within the pulsar magnetosphere itself. We can use the absence of any additional delays in our timing fits to constrain the composition of the magnetosphere and the emission height above the magnetic poles.

6.1. Aberration and retardation

According to radius-to-frequency mapping (Ruderman & Sutherland 1975; Manchester & Taylor 1977; Cordes 1978), high frequency emission is thought to originate close to the neutron star whilst low frequency emission comes from higher in the pulsar magnetosphere. The range of emission heights at different frequencies can be obtained, independently of the radius-to-frequency mapping model, through the effects of differential aberration and retardation (Cordes 1978).

Retardation is the time delay caused by the difference in path length from the different emission sites. For emission which originates at an altitude r_min < r_max, the time delay between the two paths is given by: (17)Aberration is caused by the co-rotation of the magnetosphere, which bends the radiation beam towards the direction of rotation and therefore causes pulses to arrive earlier than they would if they were to travel along a straight path. Aberration increases at larger radii, so the result is to delay pulses from low altitudes by: (18)where α is the angle between the pulsar’s magnetic and rotational axes. Combining the two effects gives the total time delay for a given radius: (19)In our data, assuming the fiducial points in our frequency-dependent models are aligned correctly, there is no remaining structure in the residuals that has not been successfully modelled by a quadratic frequency dependence. This can be used to constrain ΔR = r_max − r_min, the distance between the heights at which different frequencies are emitted (see Table 6).

For pulsars B0329+54 and B1133+16, which exhibit typical “conal” behaviour (see for example Rankin 1983a), it is also possible to use radius-to-frequency mapping to determine upper limits on the absolute height of emission from the pulsar surface. At low frequencies radius-to-frequency mapping suggests that pulse profiles are broadened as the star’s dipolar magnetic field lines move further apart high in the magnetosphere. For a neutron star with a dipolar magnetic field the ratio of the widths of the profiles (θ_ν2 > θ_ν1) at two frequencies (ν₁ > ν₂) can be related to the emission heights by the equation (Cordes 1978): (20)This equation can be used in conjunction with the values of ΔR derived earlier to determine upper limits on r_max, the absolute height of the 40 MHz emission. This analysis was not performed for pulsars B0809+74 and B1919+21, as our observations of their pulse profile evolution do not agree with the standard picture of radius-to-frequency mapping (see Sect. 7).

Our limits agree well with previous papers, such as those by Cordes (1978); Matese & Whitmire (1980); Karuppusamy et al. (2011), who also failed to detect the effects of aberration and retardation in low frequency pulsar timing data. It is also interesting to compare our findings to those of Kramer et al. (1997) who performed a similar analysis at high frequencies (1.4–32 GHz) and found that r_max < 310 km for PSR B1133+16 and r_max < 320 km for PSR B0329+54.

Our limits of r_max < 110 km for PSR B1133+16, significantly improve upon the previous low frequency (<200 MHz) limits for PSR B1133+16 set by Cordes (1978) and Matese & Whitmire (1980), who found r_max < 630 km; and Karuppusamy et al. (2011) who improved upon this, finding r_max < 560 km. This is predominantly because we have a large fractional bandwidth and high sensitivity at low frequencies. Our limits are also improved by the frequency-dependent models that were used, which allowed us to test how well our fiducial points fit the data, reducing the uncertainty in Δt_A/R significantly.

The uncertainties in our measurements are dominated by the uncertainties associated with α, which unfortunately, are not well constrained (see Everett & Weisberg 2001, for a discussion). No uncertainties on α were given in Lyne & Manchester (1988), and because of this, the uncertainties on ΔR and r_max are impossible to determine definitively. However, as α only appears in Eq. (19) through a factor of (1 + sinα) our measurements should be within a factor of two of the true value.

The implication of these limits is that pulsar emission from all radio frequencies is produced inside a very small region of the magnetosphere⁸. All of the radio emission from B1133+16 comes from within 11 stellar radii (using the canonical neutron star radius of 10 km, as used in Kramer et al. 1997), a fact which could have implications for future models of the pulsar magnetosphere.

6.2. Magnetospheric propagation effects

The pulsar magnetosphere is a complex system with strong magnetic fields and high concentrations of relativistic charged particles. As high frequency emission is thought to originate closer to the neutron star surface, it has more of the magnetosphere to travel through. This means that, under the assumption of radius-to-frequency mapping, we might expect to measure a slightly higher value for DM at high frequencies than at low frequencies, changing the way dispersion delay scales with frequency.

We see no evidence for any deviation from a ν^-2 power law in our data, suggesting that dispersion from within the magnetosphere is either not present, too small to detect, or indistinguishable from the cold plasma dispersion law (at least in the fiducial components). This suggests that either the column density of the plasma in the magnetosphere is too small to cause refraction, or emission can somehow escape the magnetosphere without being refracted. This could be because the emission of the fiducial component propagate via the extraordinary mode, which does not suffer refraction (Barnard & Arons 1986). In the pulsar magnetosphere, the electrons are very tightly bound by the magnetic field lines and cannot move transverse to their direction. This means that photons cannot excite them, effectively making them invisible and setting their refractive index to 1 (see for example Michel 1991). It should be noted that other modes of propagation, which can be refracted (ordinary modes) are also possible in the magnetosphere, although from this evidence, they do not seem to be present in the fiducial components.

6.3. Super-dispersion

Super-dispersion was proposed by Shitov & Malofeev (1985) after they observed that the DMs of low frequency pulsar observations seemed to be systematically higher than the DMs of the same pulsars at high frequencies.

They explained the delay by magnetic sweepback (Shitov 1983), in which a pulsar’s magnetic field lines get bent backwards in the opposite direction to the spin of the neutron star. This causes emission from higher up in the magnetosphere (at low frequency) to reach us slightly later than the corresponding emission from closer to the neutron star surface (high frequency). This model was supported by further evidence from Kuz’min (1986) who observed super-dispersion in eight more pulsars and observed a correlation between 1/P and the observed super-dispersive delay and also Shitov et al. (1988) who observed the effect in pulsars B0834+06, B1133+16, B1508+55, B0832+26 and B1642 − 03, and also noted that the delay appeared to be time variable in B1133+16 and B0809+74.

In a later paper (Kuz’min et al. 2008) this was cast into doubt as the super-dispersion in Crab giant pulses corresponded to more than 1 period, a delay which cannot be explained by the twisting of magnetic field lines in the pulsar magnetosphere. Super-dispersion was also not seen by Phillips & Wolszczan (1992) and Ahuja et al. (2005), who did, however see an excess delay at high frequencies.

We see no evidence for extra delays at low frequencies and can place a limit on the super-dispersive delay of  ≲ 1 ms at 40 MHz for the pulsars in our sample. We speculate that the delay which was seen in PSR B0809+74 was actually due to either pulse profile evolution or an incorrect fiducial point (see Sect. 7).

7. Pulse profile evolution

7.1. Impact on pulsar timing

We have shown that it is possible to find an analytical fiducial point for each of the pulsars in our sample, which is valid between 48 MHz and 1780 MHz, and which satisfies the ν^-2 frequency dependence expected from the dispersive delay. Building a frequency-dependent template around this fiducial point to include the effects of pulse profile evolution across the band significantly reduces the systematic errors caused by pulse profile evolution, and improves the precision of timing observations.

Pulse profile evolution can cause systematic errors in pulse arrival times which can be on the order of milliseconds for normal (slow) pulsars. The size of this effect is largest when the profile evolution is asymmetric, as noted in Ahuja et al. (2007), but it still plays a role in relatively symmetric profiles like that of PSR B0329+54. In reality there are very few pulsars that have truly symmetric pulse profile evolution, so in order to obtain ever higher timing precision (on the order of microseconds for normal pulsars or tens of nanoseconds for millisecond pulsars) it is crucial to account for the evolution of pulse shape across the band.

Although the frequency evolution across the relatively small bandwidths used up until recently will limit the effect in an individual band, it will still manifest when one tries to combine data from more than two frequencies as it cannot be absorbed into a fit for a dispersion delay. The problem becomes more acute as we search for greater sensitivity by using wider and wider instantaneous bandwidths. Here the determination of the time of arrival itself is affected by the evolution of the profile and the dispersion delay. The method presented here provides a way in which one can use very wide band data to model the profile sufficiently to build a template which incorporates all these effects, although it remains to find the optimal way to extract a time of arrival from these data.

7.2. Models

To address the profile evolution induced errors in the residuals it was necessary to make frequency-dependent models for each of the pulsars in our sample. We were able to model the profile evolution of the four pulsars over seven octaves of frequency using analytic models of Gaussian fits to the data. Where necessary (in PSR B0329+54) the models were also made to include the effects of interstellar scattering. This was modelled by convolving an exponential tail (whose length was fit to the particular frequency band) with the Gaussian components. Although the models used are simple, they describe the shape of the pulse profile very well (as shown below), and are very effective in reducing the systematic errors seen in our timing residuals as a function of frequency. Our timing residuals were used to test the validity of each of our models by determining how well the frequency-dependent templates remove the different systematic errors attributed to profile evolution in each observation. The power-law-dependencies of all of the fitted parameters in the models are given in Table 7.

In all cases, the evolution of the relative amplitudes of the peaks was very difficult to fit. In general the peak heights within an observing band could be approximated well by a power law, but the parameters of each power law varied significantly from one observing band to the next, leading to discontinuities in our model. This could be because the pulse components do not have the same spectral index across the entire frequency range (i.e. they have one or more spectral breaks), which could cause the observed discontinuities in the “ratio of peaks” parameter. The only means to test this hypothesis requires accurate flux measurements for the various components and since absolute flux calibration is not yet available to us, the present data set cannot be used to provide a conclusive answer to this particular question.

7.2.1. PSR B0329+54

Gangadhara & Gupta (2001, hereafter GG) model PSR B0329+54 as a central core component surrounded by four nested conal rings. Our model agrees well with this; we use a central core and two cones (five components) because rings 2 and 4 of the GG model are too weak for us to detect with any of our instruments. Figure 9 shows our model.

The core component (component 3), is the fiducial point for the observation and so we have set it to be constant in height and position, allowing the other components to vary. The outer cone (components 1 and 5) fades at low frequencies, and is too weak to see in the LOFAR observations. Interestingly, the two sides of the cone fade at different rates, component 5 has a steep spectrum, and is not detected in the HBA or LBA observations (it was detected by GG at 606 MHz and 325 MHz), whilst component 1 fades more gradually, and does not disappear until the LBA observations. The inner cone (components 2 and 4) is brighter and is visible all the way down to the LBAs (although the scattering tail makes it more difficult to model the components here). Again there is evidence for the two sides of the cone showing different spectral indices. Component 2 fades more than component 4 in the high frequency observations, but component 4 fades more in the LBA band.

The widths of all of the components seem to remain remarkably constant over the entire frequency range of our observations. In fact, PSR B0329+54 can be modelled with constant component widths from frequencies between 40 MHz and 8 GHz. It is difficult to say conclusively whether this model reflects a feature which is intrinsic to the pulsar, as the components all overlap, making them difficult to model. However a model with fixed component widths is simpler and this makes it much easier to track how the pulse profile evolves at different frequencies, how the components move and how their brightnesses vary in relation to each other. It is also worth noting that similar behaviour has been found in the Vela pulsar by Keith et al. (2011).

Both cones are very asymmetrical in terms of the relative brightness of their two peaks and their relative positions in relation to the central component (see Fig. 10). Compared to a model where the cones are symmetric around a central component, the outer cone is skewed by approximately 5 degrees towards earlier pulse longitudes. The leading component moves away from the central component with decreasing frequency, whilst the trailing component seems to move slightly closer. The inner cone is also skewed by approximately 5 degrees, but in the opposite direction. The components both move out from the main pulse, but component 4 moves out quite quickly, and replaces component 5 as the “postcursor” in the HBAs, whilst component 2 remains in roughly the same place until the top of the LBAs, when it begins to move out.

Our model agrees well with the GG model of a central core component surrounded by cones. The cones, however, show asymmetrical behaviour which is not expected for a dipolar magnetic field, the radio emission should come from concentric circles centred on the magnetic pole (see for example Oster & Sieber 1976).

7.2.2. PSR B0809+74

PSR B0809+74 is a rather controversial pulsar, which has shown evidence of an “absorption feature” (Bartel 1981; Rankin 1983b) and was also one of the candidates for “super-dispersion” (Shitov & Malofeev 1985, see Sect. 6.3).

Absorption was proposed by Bartel et al. (1981) who noticed that the DM found in 102 MHz observations was significantly different from the value found from observations at 1720 MHz. The reason for this was that their fiducial point was on the leading edge of the low frequency pulse profile, and the trailing edge of the high frequency pulse profile. It appeared that the low frequency pulse profile was “missing” radiation from the leading edge, which they suggested, was removed by cyclotron absorption. Further evidence for this model was also provided in a subsequent paper (Bartel 1981), which found that the profile of B0809+74 gets significantly narrower below 1 GHz. Rankin (1983b) found similar absorption features in at least eight other pulsars.

PSR B0809+74 has two overlapping components, which are normally thought to be conal. In accordance with this thinking, we fitted the data from the simultaneous observations (marked with arrows in Fig. 11) with two Gaussian components and set our fiducial point as the midpoint of the profile. This model produced large systematic errors in the TOAs at different frequencies. On a closer examination of the profiles, the reason for the timing errors became apparent, the separation of the components cannot be modelled as a simple power law. In the observations at 1400 MHz the components get closer together as frequency decreases, whereas in the low frequency data they move further apart.

In a second model, we used three components, one central component, a precursor and a postcursor. The narrower, central component is taken as the fiducial point of the profile. At high frequencies the precursor (component 1 in Fig. 11) moves towards the central component. Somewhere in the frequency range 200–1000 MHz (which was not present in the simultaneous observations) this component fades. Then, in the low frequency data, the third component appears and begins to move away from the central component, towards later pulse phase.

The precursor and the postcursor in this model have the same width and their positions can both be modelled by a single power law (see Fig. 12). This suggests that the two components may instead be a single component, which drifts through the pulse profile.

In our final model (see Fig. 11), we used two Gaussian components, a narrow component (component 2), which is the fiducial point of the pulse profile, and a broader component (component 1), which starts as a precursor in the high frequency observations and drifts through the pulse profile, arriving at a later phase (as a postcursor) at low frequencies. Using the narrow component as a fiducial point removes the systematic errors from TOAs, which provides strong verification of this model.

Further evidence in favour of this model comes from archival pulse profiles from the European Pulsar Network (EPN) database (Gould & Lyne 1998)⁹. These pulse profiles (at 410 MHz, 606 MHz and 925 MHz) are also shown in Fig. 11, along with an interpolation of our model to these frequencies. We have allowed the relative heights of the components to vary, but their positions and widths are determined by our model. Without prior knowledge about this frequency range, our model accurately predicts the shape of the pulse profile.

The fitted solution for the profile at these frequencies (where the two components overlap), is not unique. However, it is obvious from considerations of the timing residuals that the midpoint of the two components is not the fiducial point of the pulse profile and the fact that such a simple model can explain the observed profile evolution of the pulsar over such a broad frequency range is compelling.

We also compare our model to some more recent LBA observations taken using the LOFAR superterp¹⁰ (LOFAR observation ID: L30803). The data quality is significantly improved because the superterp has roughly three times the collecting area of DE601, and the delays between the dipoles have recently been calibrated. The pulse profiles are made from 6 MHz segments of bandwidth between 15 MHz and 57 MHz. Our model accurately describes the two components down to roughly 40 MHz, where the broader component begins to move away from the central component more slowly. This could be due to a mode change or some more complex pulse profile evolution at low frequencies (perhaps betraying one of the magnetospheric effects discussed above).

The polarisation of the EPN profiles also shows some evidence that one of the components in our model is linearly polarised (see Fig. 13). The first component in the data is linearly polarised at high frequencies, but as frequency decreases the linear polarisation moves towards later pulse longitudes. The early LOFAR polarisation data (obs ID: L24117), shown in the bottom panel of the figure, shows that the polarisation has moved towards the latter portion of the pulse profile at 136 MHz, and arrives after the main pulse.

One argument against this model is the step in subpulse phase which is observed at 1380 MHz (see for example Edwards & Stappers 2003a). If the precursor at 1380 MHz is the same component that appears on the trailing edge of the profile at low frequencies, it is difficult to explain what happens to the phase jump, which is not observed at 328 MHz. The mechanism by which the two components are seen to move through each other, is also still unknown. It could be due to a retardation effect which is only present in one component, or refraction from the magnetosphere, but is not immediately obvious what could be causing the profile to evolve as it appears to. Further investigation into the single pulses of PSR B0809+74 at low frequencies could help to answer these questions.

It is interesting to note that the fiducial point in our model matches the fiducial point used by Bartel et al. (1981), which led the authors to speculate that part of the low frequency pulse profile was missing. Our model, shows that the radiation is not missing, but has been displaced somehow, to later pulse longitudes. Our model also elegantly explains the narrowing of the pulse profile (also attributed to absorption) discussed in Bartel (1981) and Rankin (1983b). The cumulative pulse starts at high frequencies as two fairly distinct components, the components get closer together as frequency decreases, reducing the apparent width of the profile. At around 400 MHz, the profiles are exactly on top of each other, and the pulse width is at a minimum. Below this frequency, the broader pulse continues to move towards later pulse phase and the profile width begins to increase again, reproducing the shape of the absorption.

Super-dispersion can also be explained by this model, if the centre of the two components was used as a fiducial point for the pulse profile, instead of the core component, as in our model then the pulse would look like it arrived later than expected at low frequencies.

7.2.3. PSR B1133+16

PSR B1133+16 is one of the prototypical examples of a “well resolved conal double” profile (Rankin 1983a). It is one of the brightest pulsars in the Northern sky, so it has been widely studied and is often used as evidence in favour of radius-to-frequency mapping (see for example Thorsett 1991). The separation between its components shows a continuous increase with decreasing frequency, which is thought to trace the dipolar shape of the pulsar’s magnetic field.

This is exactly what we see in our model. It has two strong components (components 1 and 3) separated by bridge emission (component 2) (see Fig. 14). The fiducial point is the midpoint between components 1 and 3, which has been defined to be the peak of the bridge emission in this model. The conal components move further apart at lower frequencies, scaling with a power law  ~ν^-0.62. This is consistent with the exponent found in Thorsett (1991) and Xilouris et al. (1996), who found exponents of  −0.50 and  −0.71 respectively. The exponent is however, significantly lower than the value in Karuppusamy et al. (2011) and Cordes (1978), who both found a power law  ~ν^-0.3.

This is because their power law fits did not include the constant term, Δθ_min, the width of the pulse profile at the surface of the neutron star. This term was proposed by Thorsett (1991), who found that the separation of the components tends towards a constant value at high frequencies. If we do not include this term in our fits, we find a power law of ν^-0.2, which is roughly consistent with those of Karuppusamy et al. (2011) and Cordes (1978). For our analysis, we used the form of the power law from Thorsett (1991), which provides a better fit to the data.

The width of the conal components as a function of frequency has been a subject of interest in the past. Mitra & Rankin (2002) found that the component widths remain constant between 40 and 3000 MHz. We also see evidence of this in our model above 80 MHz, although the component widths begin to broaden below this value. Mitra & Rankin (2002) also noticed this broadening, and attributed it to dispersive smearing across a frequency channel or scattering from the ISM. In our observations, the dispersive smearing at 48 MHz (across a single 12 kHz channel) is  ~1.5^o, which is enough to explain the observed broadening in the profile¹¹. However, even disregarding this low frequency broadening Mitra & Rankin (2002) found that the spacing between the components between 100 MHz and 10 GHz changes too rapidly to be caused by a dipolar magnetic field.

PSR B1133+16 is the pulsar which is most consistent with radius-to-frequency mapping of all the pulsars in our sample. However, in Sect. 6 we showed that its emission is confined to a very narrow region in the magnetosphere (<59 km) which is incompatible with the standard radius-to-frequency model. Radius-to-frequency mapping assumes that the emission at a given emission height traces the last open field line in the pulsar magnetosphere. From geometrical arguments (see for example Lorimer & Kramer 2005) it is possible to write the opening angle of the last open field line of a non-relativistic dipole, ρ, as a function of emission height, r: (21)where R_LC is the radius of the light cylinder of the pulsar. ρ is the maximum value that it is possible for components to be separated by as each component originates from one side of the dipolar field. The maximum increase in component separation predicted by this equation for PSR B1133+16 over the height range, ΔR, derived earlier is 2.8°. In the pulse profile of PSR B1133+16 the components are seen to move apart by 6.3    ±    0.5° between 1780 MHz and 48 MHz. This, coupled with the further evidence by Mitra & Rankin (2002) suggests that PSR B1133+16 cannot be explained through radius-to-frequency mapping in a simple dipolar magnetic field.

7.2.4. PSR B1919+21

PSR B1919+21 was the first pulsar ever discovered (Hewish et al. 1968) and is often referred to as a classic radio pulsar. In fact, PSR B1919+21 disagrees with the classic picture of a pulsar in almost all aspects of its pulse profile evolution.

Lyne et al. (1971) showed that the two components of the pulsar get closer together as frequency decreases in the range 150–3000 MHz, which is the opposite to what one would expect from radius-to-frequency mapping. This result was confirmed by Sieber et al. (1975) who also showed that the pulse profile seems to get broader again below 150 MHz. Mitra & Rankin (2002) found that the width of the profile and separation of the components was approximately constant in their observations between 50 and 5000 MHz. They suggested that the lack of any radius-to-frequency mapping could be due to the emission originating from an “inner cone”, which is located closer to the “core” emission component, and so the effects of radius-to-frequency mapping are less pronounced.

In our model (see Fig. 3) the two components move closer together at low frequencies (in agreement with Lyne et al. 1971). We did not, however, see the component separation increasing below 150 MHz (as reported by Sieber et al. 1975). This behaviour is similar to that of PSR B0809+74 at high frequencies, although in this case the components are only expected to pass through each other at  ~1 MHz. Another curious feature of our model is that whilst component 2 gets broader at low frequencies (as expected), component 1 appears to get narrower. Radius-to-frequency mapping suggests that both components should get broader at low frequencies, as the emission region gets wider. In an attempt to conform with this idea, we tried to fit a three-component model. This fit the pulse profile reasonably well, and also agreed better with the radius-to-frequency model; the components didn’t move closer together, and stayed at a constant width. However, a suitable fiducial point could not be determined, and there were large systematic errors in our timing residuals so the model was rejected.

8. Discussion

8.1. Profile evolution

We have shown that pulse profile evolution can introduce large errors into pulsar timing data, in agreement with the work of Ahuja et al. (2007). These errors can be as large as a few milliseconds in some cases, depending on the period of the pulsar, and how asymmetric the pulse profile is. A frequency-dependent model of the pulse profile can be used to reduce these errors. Using this method, it was possible to define an analytical fiducial point in the pulse profile of each of the pulsars in our sample. This fiducial point is valid to within a few milliseconds (corresponding to  ~1 degree in pulse phase), although the model does not remove the timing errors completely. Small timing errors remain because the model is not an exact fit to the observed profile, and subtle differences between the shape of the modelled templates and the data lead to systematics in the cross-correlation. Further work needs to be done to explore how to better remove these timing errors, and how they vary with time.

We have found that radio emission from all of the pulsars in our sample originates from a narrow range of heights in the magnetosphere. The narrow ranges found do not fit well with models of radius-to-frequency mapping in a dipolar magnetic field. In addition, all of the pulsars in our sample show at least one other trait which cannot be explained by radius-to-frequency mapping.

The asymmetric cones which we observe in PSR B0329+54 were also observed by GG, who attributed their asymmetry to aberration and retardation. However, we have not detected any aberration or retardation effects in our timing residuals and we also find that the emission from the inner cone seems to be concentrated to within 183 km of the neutron star surface. The (~5 degree) skew in the cone corresponds to a time difference of  ~10 ms, which is much greater than the aberration and retardation effects possible from within this height range, which are  ~1 ms (from Eq. (19)). The fact that the outer cone is skewed in the opposite direction to the inner cone, also suggests that this cannot be explained by the standard model of a pulsar.

PSR B0809+74 has a component which starts out as a precursor at high frequencies and then drifts through the centre of the pulse profile, swapping sides with the central component and appearing as a postcursor at low frequencies. The frequency dependence of the position of the drifting component suggests that either refraction or some relationship between frequency and height (a change in height could explain a component being delayed) significantly influences one component, but is not seen in the other.

PSR B1133+16 shows emission from a very narrow range of heights and, as Mitra & Rankin (2002) showed, component separation which increases too rapidly to be produced by dipolar field lines. One explanation for this could be that there are other mechanisms at work, which act together with the traditional picture of a pulsar, complicating pulse profile evolution.

PSR B1919+21 has a profile whose width decreases at lower frequencies. This is the exact opposite of what is predicted by radius-to-frequency mapping, and so is very difficult to explain using the standard picture of a pulsar. Again, there is a clear relationship between pulse shape and frequency, but it does not seem to be explainable by radius-to-frequency mapping.

The fact that none of the pulsars in our sample behave as predicted by radius-to-frequency mapping suggests that a more complicated model of the pulsar magnetosphere is needed to describe pulse profile evolution. Although radius-to-frequency mapping has been successful in explaining some of the features seen in pulse profiles, it is clear that it cannot be used to fully describe any of the pulsars in our sample. There are, however, alternative theories which could potentially provide good fits to observational data.

In the family of models proposed by Barnard & Arons (1986), and developed further by Petrova (2000), Weltevrede et al. (2003) and Beskin & Philippov (2011), the frequency-dependent profile evolution seen in pulsars is explained by propagation effects in the pulsar magnetosphere. In these models, the radiation originates from a small region in the magnetosphere, and refraction, dispersion and the different propagation modes (i.e. extraordinary and ordinary) in the magnetosphere are responsible for the frequency evolution of the different components which are observed in pulse profiles.

Karastergiou & Johnston (2007) also provide an interesting empirical model of the pulsar magnetosphere, which could be used to explain all of the features which we have observed. They postulate that all radio emission originates from a patchy cone bounded by the last open field lines and that emission can come from any height, independently of frequency. Complex pulse profiles can then be explained by invoking emission from a range of heights, rather than assuming that the pulse profile probes the longitudinal shape of the beam at one single height.

What causes pulse profile evolution is still an important question, and will be vital in understanding the pulsar emission mechanism, and for studies of pulsar geometries in the future. At this stage, it is still difficult to discriminate between the many models that exist, but next generation telescopes, like LOFAR, will be excellent tools for studying this effect.

8.2. Magnetospheric effects

The argument that a more sophisticated model is needed to describe radio emission from pulsars is also supported by considerations of aberration and retardation effects on our data. From these arguments, we have shown that radio emission from all of the pulsars in our sample is confined to a very small region in the pulsar magnetosphere, which supports the ideas of Barnard & Arons (1986) and Petrova (2000).

We have also shown that, as there is no departure from a ν^-2 dispersion law in our data, there is no evidence for super-dispersive delays or refraction from within the pulsar magnetosphere in any of the pulsars in our sample. However, whilst we don’t see a frequency-dependent delay in the timing residuals, refraction may be needed to explain the broad component of PSR B0809+74, which is seen to drift through the pulse profile.

8.3. ISM effects

In our data, we see no evidence for any deviation from the cold plasma dispersion law, suggesting that second-order ISM effects in these pulsars introduce additional time delays  ≲ 50 ns at normal pulsar timing frequencies (1400 MHz). The parameters of the pulsars in our sample are typical of those found in the PTAs, the only difference being their longer pulse periods. This suggests that the ISM may not cause as much of a hindrance to pulsar timing projects as first feared (Hemberger & Stinebring 2008).

The fact that no unexpected delay was detected in any of our observations also means that (at least along these lines of sight) the ISM appears to be relatively smooth, with no large, dense structures. These findings agree with the idea that scattering is dominated by one or two small, but relatively high density regions as discussed by Stinebring (2006) and Brisken et al. (2010). We have determined an upper limit on , which can be used as an indicator as to how much the DM is likely to change in the future. Comparing these predictions with reality will be a useful cross-check of how well this relation works. We have also been able to place upper limits on the scattering measure and the emission measure. Although these limits are weak compared to NE2001 (Cordes & Lazio 2002), pulsar timing could provide an independent method of measuring both of these quantities in the future.

8.4. Future observations

The constraints set in this paper will be improved significantly by taking similar observations when LOFAR is completed, using more stations at a lower observing frequency. The sensitivity of LOFAR has already improved by a factor of five since the observations for this paper were taken, and it is expected that it will increase significantly again soon, when the core stations can all be combined coherently. Increased sensitivity, particularly in the LBAs, will reduce the error bars seen at low frequencies in our timing residuals, which dominate the uncertainty in our measurement.

We could also increase our precision by observing at lower frequencies. LOFAR will soon be able to routinely observe with high sensitivity at frequencies as low as 15 MHz (see Fig. 11), where the second-order ISM delays are expected to be at least an order of magnitude larger than they are at 40 MHz.

By repeating this experiment in the future on the same set of pulsars, we could test whether pulse profile evolution is stable with time, and also track variations in the DM with great accuracy. Both of these parameters are not completely understood, and are vital for high-precision pulsar timing. It would also be interesting to perform this experiment on a millisecond pulsar. A faster rotation rate and a narrower pulse (in absolute terms) means that TOAs can be determined more accurately, which would improve our constraints by at least an order of magnitude.

---

Acknowledgments

We would like to thank Jim Cordes for his insight and useful discussions, Christine Jordan for arranging the observations from Jodrell Bank, and the anonymous referee for their insightful comments. LOFAR, the LOw Frequency ARray designed and constructed by ASTRON, has facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the International LOFAR Telescope (ILT) foundation under a joint scientific policy. This publication made use of observations taken with the 100-m telescopes of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg. Ben Stappers, Patrick Weltevrede and the Lovell observations are supported through an STFC rolling grant. Tom Hassall is the recipient of an STFC studentship. Jason Hessels is a Veni Fellow of the Netherlands Foundation for Scientific Research. Joeri van Leeuwen and Thijs Coenen are supported by the Netherlands Research School for Astronomy (Grant NOVA3-NW3-2.3.1) and by the European Commission (Grant FP7-PEOPLE-2007-4-3-IRG #224838). Aris Karastergiou is grateful to the Leverhulme Trust for financial support. Joris Verbiest is supported by the European Union under Marie Curie Intra-European Fellowship 236394. Charlotte Sobey is supported by the DFG (German Research Foundation) within the framework of the Research Unit FOR 1254, Magnetisation of Interstellar and Intergalactic Media: The Prospects of Low-Frequency Radio Observations.

References

All Tables

All Figures