© ESO, 2016

1. Introduction

The cumulative (i.e. average) pulse profiles of radio pulsars are the sum of hundreds to thousands of individual pulses, and are, loosely speaking, a unique signature of each pulsar (Lorimer 2008). They are normally stable in their morphology and are reproducible, although several types of variation have been observed both for non-recycled pulsars (see Helfand et al. 1975; Weisberg et al. 1989; Rathnasree & Rankin 1995 and Lyne et al. 2010) and for millisecond pulsars (Liu et al. 2012). For most pulsars, this cumulative pulse profile morphology often varies (sometimes drastically, sometimes subtly) as a function of observing frequency because of a number of intrinsic effects (e.g. emission location in the pulsar magnetosphere) and extrinsic effects (i.e. due to propagation in the interstellar medium; ISM), see for instance Cordes (1978). As Fig. 1 shows, pulse profile evolution can become increasingly evident at the lowest observing frequencies (<200 MHz).

Mapping profile evolution over a wide range of frequencies can aid in modelling the pulsar radio emission mechanism itself (see e.g. Rankin 1993 and related papers of the series) and constraining properties of the ISM (e.g. Hassall et al. 2012). Since many of the processes that affect the pulse shape strongly depend on the observing frequency, observations at low frequencies provide valuable insights on them. The LOw-Frequency ARray (LOFAR) is the first telescope capable of observing nearly the entire radio spectrum in the 10−240 MHz frequency range, the lowest 4.5 octaves of the “radio window” (van Haarlem et al. 2013) and (Stappers et al. 2011).

Low-frequency pulsar observations have previously been conducted by a number of telescopes, (e.g. Gauribidanur Radio Telescope (GEETEE): Asgekar & Deshpande 1999; Large Phased Array Radio Telescope, Puschino: Kuzmin et al. (1998), Malov & Malofeev 2010; Arecibo: Hankins & Rankin 2010; Ukrainian T-shaped Radio telescope, second modification (UTR-2): Zakharenko et al. 2013), and simultaneous efforts are being undertaken by other groups in parallel (e.g. Long Wavelength Array (LWA): Stovall et al. 2015, Murchison Widefield Array (MWA): Tremblay et al. 2015). Nevertheless, LOFAR offers several advantages over the previous studies. Firstly, the large bandwidth that can be recorded at any given time (48 MHz in 16-bit mode and 96 MHz in 8-bit mode) allows for continuous wide-band studies of the pulse profile evolution, compared to studies using a number of widely separated, narrow bands (e.g. the 5 × 32 and 20 × 32 kHz bands used at 102 MHz in Kuzmin et al. 1998). Secondly, LOFAR’s ability to track sources is also an advantage, as many pulses can be collected in a single observing session instead of having to combine several short observations in the case of transit instruments (e.g. Izvekova et al. 1989). This eliminates systematic errors in the profile that are due to imprecise alignment of the data from several observing sessions. Thirdly, LOFAR can achieve greater sensitivity by coherently adding the signals received by individual stations, giving a collecting area equivalent to the sum of the collecting area of all stations (up to 57 000 m² at HBAs and 75 200 m² at LBAs, see Stappers et al. 2011). Finally, LOFAR offers excellent frequency and time resolution. This is necessary for dedispersing the data to resolve narrow features in the profile. LOFAR is also capable of coherently dedispersing the data, although that mode was not employed here.

As mentioned above, there are two types of effects that LOFAR will allow us to study with great precision. One of these are extrinsic effects related to the ISM. Specifically, the ISM causes scattering and dispersion. Mean scatter-broadening (assuming a Kolmogorov distribution of the turbulence in the ISM) scales with observing frequency as ν^-4.4. The scattering causes delays in the arrival time of the emission at Earth, which can be modelled as having an exponentially decreasing probability of being scattered back into our line-of-sight: this means that the intensity of the pulse is spread in an exponentially decreasing tail. Dispersion scales as ν^-2 and is mostly corrected for by channelising and dedispersing the data (see e.g. Lorimer & Kramer 2004). Nonetheless, for filterbank (channelised) data some residual dispersive smearing persists within each channel: ${\mathit{t}}_{\mathrm{DM}}\mathrm{=}\mathrm{8.3}\mathrm{·}\mathrm{DM}\frac{\mathrm{\Delta }\mathit{\nu }}{{\mathit{\nu }}^{\mathrm{3}}}\mathit{\mu }\mathrm{s}\mathit{,}$(1)where DM is the dispersion measure in cm^-3 pc, Δν is the channel width in MHz, and ν the central observing frequency in GHz. This does not significantly affect the profiles that we present here (at least not at frequencies above ~50 MHz) because the pulsars studied have low DMs and the data are chanellised in narrow frequency channels (see Sect. 2). Second-order effects in the ISM may also be present, but have yet to be confirmed. For instance, previous claims of “super-dispersion”, meaning a deviation from the ν^-2 scaling law (see e.g. Kuz’min et al. 2008 and references therein), were not observed by Hassall et al. (2012), with an upper limit of ≲50 ns at a reference frequency of 1400 MHz.

The second type of effects under investigation are those intrinsic to the pulsar. One of the most well-known intrinsic effects are pulse broadening at low frequencies, which has been observed in many pulsars (e.g. Hankins & Rickett 1986 and Mitra & Rankin 2002), while others show no evidence of this (e.g. Hassall et al. 2012). One of the theories explaining this effect is radius-to-frequency mapping (RFM, Cordes 1978): it postulates that the origin of the radio emission in the pulsar’s magnetosphere increases in altitude above the magnetic poles towards lower frequencies. RFM predicts that the pulse profile will increase in width towards lower observing frequency, since emission will be directed tangentially to the diverging magnetic field lines of the magnetosphere that corotates with the pulsar. An alternative interpretation (McKinnon 1997) proposes birefringence of the plasma above the polar caps as the cause for broadening: the two propagation modes split at low frequencies, causing the broadening, while they stay closer together at high frequencies, causing depolarisation (this is investigated in the LOFAR work on pulsar polarisation, see Noutsos et al. 2015).

In this paper, we present the average pulse profiles of 100 pulsars observed in two LOFAR frequency ranges: high band (119–167 MHz, 100 profiles) and low band (15–62 MHz, 26 out of the 100 profiles). We compare the pulse profile morphologies with those obtained around 350 and 1400 MHz with the WSRT and Lovell telescopes, respectively, to study their evolution with respect to a magnetospheric origin and DM-induced variations. We do not discuss here profile evolution due to the effects of scattering in the ISM, which will be the target of a future work (Zagkouris et al., in prep.). In Sect. 2 we describe the LOFAR observational setup and parameters. In Sect. 3 we describe the analysis. In Sect. 4 we discuss the results, and in Sect. 5 we conclude with some discussion on future extensions of this work.

2. Observations

The observed sample of pulsars was loosely based on a selection of the brightest objects in the LOFAR-visible sky (declination >−30°), using the ATNF Pulsar Catalog¹ (Manchester et al. 2005) for guidance. Because pulsar flux and spectral indices are typically measured at higher frequencies, we also based our selection on the previously published data at low frequencies (Malov & Malofeev 2010). Since the LOFAR dipoles have a sensitivity that decreases as a function of the zenith angle, all sources were observed as close to transit as possible, and only the very brightest sources south of the celestial equator were observed.

We observed 100 pulsars using the high-band antennas (HBAs) in the six central “Superterp” stations (CS002−CS007) of the LOFAR core². The observations were performed in tied-array mode, that is, a coherent sum of the station signals using appropriate geometrical and instrumental phase and time delays (see Stappers et al. 2011 for a detailed description of LOFAR’s pulsar observing modes and van Haarlem et al. 2013 for a general description of LOFAR). The 119−167-MHz frequency range was observed using 240 subbands of 195 kHz each, synthesised at the station level, where the individual HBA tiles were combined to form station beams. Using the LOFAR Blue Gene/P correlator, each subband was further channelised into 16 channels, formed into a tied-array beam. The linear polarisations were summed in quadrature (pseudo-Stokes I), and the signal intensity was written out as 245.76 μs samples. The integration time of each observation was at least 1020 s. This was chosen to provide an adequate number of individual pulses, so as to average out the absolute scale of the variance associated with pulse-phase “jitter” to the cumulative profile. The jitter, also termed stochastic wide-band impulse modulated self-noise (SWIMS, as in Osłowski et al. 2011), is the variation in individual pulse intensity and position with respect to the average pulse profile (see also Cordes 1993 and Liu et al. 2012 and references therein). This variation does not significantly affect the measurements that have been carried out for the scope of this paper (i.e. pulse widths, peak heights), but we have checked that the resulting profile was stable on the considered time scales by dividing each observation into shorter sections and comparing the shapes of the resulting profiles with the overall profile. In the cases where stability was not achieved, we used longer integration times. Regardless, in almost all cases the profile evolution with observing frequency is a significantly stronger effect at low frequencies.

Twenty-six of the brightest pulsars were also observed using the Superterp low-band antennas (LBAs) in the frequency range 15−62 MHz. To mitigate the larger dispersive smearing of the profile in this band, 32 channels were synthesised for each of the 240 subbands. The sampling time was 491.52 μs. The integration time of these observations was increased to at least 2220 s to somewhat compensate for the lower sensitivity at this frequency band (e.g. because of the higher sky temperature).

For some sources, 17-minute HBA observations with the Superterp were insufficient to achieve acceptable signal-to-noise (S/N) profiles. For these, longer integration times (or more stations) were needed. Hence, some of the pulsars presented here were later observed with 1 hr pointings as part of the LOFAR Tied-Array All-Sky Survey for pulsars and fast transients (LOTAAS³: see also Coenen et al. 2014), which commenced after the official commissioning period, during Cycle 0 of LOFAR scientific operations, and is currently ongoing. LOTAAS combines multiple tied-array beams (219 total) per pointing to observe both a survey grid as well as known pulsars within the primary field-of-view.

3. Analysis

The LOFAR HDF5⁴ (Hierarchical Data Format5, see e.g., Alexov et al. 2010) data were converted to PSRFITS format (Hotan et al. 2004) for further processing. Radio frequency interference (RFI) was excised by removing affected time intervals and frequency channels, using the tool rfifind from the PRESTO⁵ software suite (Ransom 2001). The data were dedispersed and folded using PRESTO and, in a first iteration, a rotational ephemeris from the ATNF pulsar catalogue, using the automatic LOFAR pulsar pipeline “PulP”. The number of bins across each profile was chosen so that each bin corresponds to approximately 1.5 ms.

The profiles obtained with the HBA and, where available, LBA bands were compared with the profiles obtained with the WSRT at ~350 MHz (from here onwards “P-band”) and at ~1400 MHz, or with the Lovell Telescope at the Jodrell Bank Observatory at ~1500 MHz (from here onwards “L-band”). The WSRT observations that we used were performed mostly between 2003 and 2004 (see Weltevrede et al. 2006, “WES” from here onwards, and Weltevrede et al. 2007 for details). The Lovell observations were all contemporary to LOFAR observations, therefore in the cases where both sets of observations were available, we chose the Lovell ones because they are closer to or overlap the epoch of the LOFAR observations. At L-band we used Lovell observations for all but three pulsars: B0136+57, B0450−18, and B0525+21. In a handful of cases, where no profile at P-band was available from WES, we used the data from the European Pulsar Network (EPN) database⁶.

To attempt to align the data absolutely, we generated ephemerides that spanned the epochs of the observations that were used. This did not include those from the EPN database, however. Ephemerides were generated from the regular monitoring observations made with the Lovell Telescope. The times of arrival were generated using data from an analogue filterbank (AFB) up until January 2009 and a digital filterbank (DFB) since then, with a typical observing cadence between 10 and 21 days. The observing bandwidth was 64 MHz at a central frequency of 1402 MHz and approximately 380 MHz at a central frequency of 1520 MHz for the AFB and DFB, respectively. The ephemerides were generated using a combination of PSRTIME⁷ and TEMPO, and in the case of those pulsars demonstrating a high degree of timing noise, up to five spin-frequency derivatives were fit to ensure white residuals and thus good phase alignment.

The L-band profiles were generated from DFB observations by forming the sum of up to a dozen observations, aligned using the same ephemerides used to align the multi-frequency data. We re-folded both the LOFAR and the high frequency data sets using this ephemeris. In general, where Lovell data were available, the ephemeris was created using about 100 days’ worth of data. For the WSRT observations, an ephemeris was created spanning, in some cases, ten years of data and ending at the time of the LOFAR observations. The timing procedure was the same as for the shorter data spans, except that astrometric parameters were fit and typically more spin-frequency derivatives were required. The epoch of the WSRT observations is specified in Table B.1 in the Notes column. Given the method we used to align the profiles, the timing solution is less accurate over these longer time spans than those constructed to align the Lovell data, but this too represents a good model, with a standard deviation (rms) of the timing residuals ≲1 ms. We aligned the profiles in absolute phase by calculating the phase shift between the reference epoch of the observations and the reference epoch of the ephemeris and applying this phase shift to each data set.

Some of the pulsar parameters derived from these ephemerides are presented in Table B.1. The first column lists the observed pulsars, the second and third columns list the spin period and period derivative of each pulsar, the fourth column is the reference epoch of the rotational ephemeris that was used to fold the data, and Cols. 5 and 6 list the epochs of the LOFAR HBA and LBA observations. In Cols. 7 and 8 two measurements for the DM are given: the first as originally used to dedisperse the observations at higher frequencies, and the second as the best DM obtained from the fit of the HBA LOFAR observations using PRESTO’s prepfold (Ransom 2001). The next three columns provide the pulsar’s spin-down age, magnetic field strength, and spin-down luminosity as derived from the rotational parameters according to standard approximations (see e.g. Lorimer & Kramer 2004):

$\begin{array}{ccc}& & \mathit{\tau }\mathrm{\left[}\mathrm{s}\mathrm{\right]}\mathrm{=}\mathrm{0.5}\mathit{P}\mathit{/}\mathit{Ṗ}\mathit{,}\\ & & \mathit{B}\mathrm{\left[}\mathrm{G}\mathrm{\right]}\mathrm{=}\mathrm{3.219}\mathrm{\times }{\mathrm{10}}^{\mathrm{19}}\sqrt{\mathit{P}\mathit{Ṗ}}\mathit{,}\\ & & \mathit{Ė}\mathrm{\left[}\mathrm{erg}\mathit{/}\mathrm{s}\mathrm{\right]}\mathrm{=}\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}\mathrm{\times }{\mathrm{10}}^{\mathrm{45}}\mathit{Ṗ}\mathit{/}{\mathit{P}}^{\mathrm{3}}\mathit{,}\end{array}$where P is measured in s and Ṗ in s s^-1. The resulting aligned profiles for the 100 pulsars are shown in Fig. B.1 in Appendix B. The star next to the name of the band (P-band in most cases, with the exception of B0136+57, where we used the P-band for absolute alignment) indicates that the corresponding band was aligned by eye based on the absolute alignment between the LOFAR data and the other high-frequency band. The alignment was made based only on the choice of a specific point along the rotational phase of the pulsar, at the reference epoch of the ephemeris, but unmodelled DM variations can be responsible for extra, albeit small, phase shifts (up to a few percent, see Table 1). Indeed, observations performed at different times, quite far apart, and at different frequencies, can possess quite different apparent DMs (up to some tenth of a percent, see Table 1). DMs that are due to the ISM have, as expected, a time dependence (see You et al. 2007; Keith et al. 2013), and these differences can become quite relevant especially at the lowest LOFAR frequencies. We chose to re-dedisperse all the profiles, LOFAR and high-frequency ones, using the DM obtained as the best DM with prepfold for the HBA LOFAR observations. prepfold determines an optimum DM by sliding frequency subbands with respect to each other to maximise the S/N of the cumulative profile. The intra-channel smearing caused by DM over the bandwidth at the centre frequency is indicated by the filled rectangle next to each profile in Fig. B.1.

Only in a few cases, documented in Table 1, was the DM obtained from the LOFAR observations not used for the alignment. Those are the cases for which, also evident from Fig. B.1, the intra-channel smearing caused by DM is a substantial fraction of the profile width (similar to or higher than the on-pulse region), and therefore the quality of the measurement is lower than that obtained at a higher observing frequency. On the other hand, in some other cases (although rare, see Table B.1), even the LBA measurement was good enough to provide a DM measurement, and in these cases we were able to use that for the alignment. In this way we obtain the best alignments, in general, even though some residual offsets could still be observed in a handful of cases.

For those cases (listed in Table 1) where the remaining offset was noticeable by eye, we investigated the possible causes after refolding and applying the new DM. We checked whether the pulsars in our sample had undergone any glitch activity during the time spanned by our ephemerides. Sixteen out of our 100 pulsars have shown glitch activity at some time. Seven of them have experienced glitches relatively close to the epoch of our observations, but only two of them during the time spanned by our ephemerides. The relevant glitch epochs of these pulsars are presented in Table 2 and were taken from the Jodrell Bank glitch archive⁸ (Espinoza et al. 2011, 2012), integrated with the ATNF pulsar archive⁹. We note that while the glitch activity could have had an influence on the shift of PSR B0525+21, the recurrent activity of PSR B0355+54 did not cause as notable an impact on the alignment. In some cases the profile is scattered in the LOFAR HBA band and rapidly becomes more scattered towards lower frequencies. This will affect the accuracy of the DM measurement, potentially causing an extra profile shift between bands. In some cases the ephemeris had a large rms timing residual, or it included higher order spin-frequency derivatives beyond the second, which is typical of unstable, “noisy” pulsars. All these cases are indicated in Table 1. In the other cases, we calculated that a DM difference ≲0.05 cm^-3 pc, compatible with the measurement uncertainties, would compensate for that shift. For this reason, we applied an extra shift to align these pulsars’ profiles by eye in Fig. B.1 and referenced this shift in Table 1.

Only a small number of pulsars in our sample show interpulses: B0823+26, B0950+08, B1822–09, B1929+10, and B2022+50. For these pulsars the profile longitude is shown in its entire phase range instead of in the phase interval 0.25−0.75, as was chosen for the other profiles. The phase-aligned profiles, which normally have their reference point at 0.5, have been shifted by + 0.25 in these cases, to help the visualisation of the interpulse. A zoom-in of the interpulse region is shown in Fig. B.2 (LBA data were excluded because none of the interpulses were detected in that band). Our sample also contains a few moding pulsars (notably B0823+26 and B0943+10). For these pulsars we caution that the profile reflects only the mode observed in the particular observations presented here. A more detailed treatment of the low-frequency profiles of B0823+26 and B0943+10 can be found in Sobey et al. (2015) and Bilous et al. (2014), respectively. In yet other cases, for instance B1857–26, no moding behaviour has previously been documented, but notable changes in the profile are seen from high to low frequencies. These might reflect different modes of emission when the various observations were taken. A more detailed discussion of some cases of peculiar profile evolution can be found in Sect. 5.3.

We fit the multi-band profiles of each pulsar using Gaussian components (see an example in Fig. 2), which are a good representation for the profiles of slow pulsars (see Kramer et al. 1994 and references therein). We used the program pygaussfit.py, of the PRESTO suite, which has the advantage of providing an interactive basis for the input parameters to the Gaussian fits. This program can be used to apply the same method as in Kramer et al. (1994; see e.g. Fig. 3 of their paper), as it provides post-fit residuals (for the discrepancy between the model and the data, see bottom half-plots of Fig. 2) that we required to have approximately the same distribution in the on-pulse as in the off-pulse region. The full rotational period, and not only the on-pulse region, was taken into account by the fit, also allowing for distinguishing interpulses or small peaks at different phases from the noise. The Gaussian components derived using this method were chosen to satisfy the condition of best fit with minimal redundancy, and no physical significance should be attributed to them. Hassall et al. (2012) have shown that it is possible to accurately model the evolution of the profile with frequency using Gaussian components, but a specific model has to be applied individually to each pulsar, requiring much careful consideration. Such an analysis is beyond the scope of this paper. In the absence of a more comprehensive treatment of scattering for LOFAR profiles, which is envisaged for subsequent papers, the most evident scattering tails of the low-frequency profiles have not been modelled and the corresponding profiles components were not included in the table and are not used in any further analysis.

We used the mathematical description of the profiles in terms of the Gaussian components to calculate the widths and amplitudes of the observed peaks. For each profile we obtained the full width at half maximum (w₅₀) and the full width at 10% of the maximum, w₁₀. To calculate the errors on the widths, we simulated 1000 realisations of each profile, using the noiseless Gaussian-based template and adding noise with a standard deviation equal to that measured in the off-pulse region of the observed profile. By measuring the widths in these realisations, we obtained a distribution of allowed widths from which we determined their standard deviation. We note that these errors are statistical only and do not take into account the validity of our assumptions, that is, the reliability of the template used. To cross-check the width of the full profile, we tried different methods. An example of how these widths were calculated is described in Appendix A and is shown in Fig. A.1. In Table B.2, Cols. 2 and 3 show w₅₀ and w₁₀ in degrees at all frequencies, calculated using the Gaussian profiles and cross-checked using the on-pulse visual inspection. The last two columns represent the calculated spectral index δ of the evolution of these widths (w₅₀ and w₁₀) with frequency as w ∝ ν^δ. The table shows that this fit can be highly uncertain. We also note that for more data, the linear fit is not always the best representation of the real trend of the profile evolution (see Fig. 3 and a more detailed discussion in Sect. 5.1.1). Column 4 of Table B.2 lists the duty cycle of each pulsar, calculated as w₁₀/P, where P is the pulsar period.

4. Results

Here we present the results of LOFAR observations of 100 pulsars, considering their profile evolution with frequency and comparing them with observations at higher frequencies. In particular, we study how the number of profile peaks, their widths, and the relative pulse phases vary with frequency.

4.1. Pulse widths

We calculated the evolution of the profile width across the frequency range covered by our observations. We chose w₁₀ for our calculations and cross-checked using w_op, as it is better suited for multi-peaked profiles than w_eff and less affected by low S/N than w_pow (see Sect. A, Fig. A.1 and Table B.2). When the measurements disagreed, the Gaussian fit was refined after visual inspection of the obtained width.

We calculated the dependence of the width of the profiles on the pulse period, considering the different frequencies separately. We note that the pulse width is not a direct reflection of the beam size or diameter (i.e. 2ρ, where ρ is beam radius). For a visual representation of the geometry see for instance Maciesiak et al. (2011) and Bilous et al. (2014). In fact, only if the observer’s line of sight cuts the emission centrally for magnetic inclination angles, α, that are not too small (i.e. α> ~ 60°), w ≈ 2ρ. In such a case, when the emission beam is confined by dipolar open field lines, we would expect a P^{− 1 / 2} dependence, which has indeed been observed when correcting for geometrical effects by transforming the pulse width into a beam radius measurement (see Rankin 1993; Gil & Krawczyk 1996; Maciesiak et al. 2012). For circular beams, profile width and beam radius are related by the relation first derived by Gil et al. (1984): ${\mathit{\rho }}_{\mathrm{10}}\mathrm{=}\mathrm{2}{\sin }^{-1}{\left[\sin \mathit{\alpha }\sin \mathrm{(}\mathit{\alpha }\mathrm{+}\mathit{\beta }\mathrm{)}{\sin }^{\mathrm{2}}\left(\frac{{\mathit{w}}_{\mathrm{10}}}{\mathrm{4}}\right)\mathrm{+}{\sin }^{\mathrm{2}}\left(\frac{\mathit{\beta }}{\mathrm{2}}\right)\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}$(5)The angle β is the impact angle, measured at the fiducial phase, φ, which describes the closest approach of our line of sight to the magnetic axis. This equation is derived under the assumption that the beam is symmetric relative to the fiducial phase. Typically, widths are measured at a certain intensity level (e.g. 50% or 10%, as here), and ρ values are derived accordingly. In many cases, profiles are indeed often asymmetric relative to the chosen midpoint, or become so as they evolve with frequency. We note that for a central cut of the beam (β = 0) and for an orthogonal rotator (α ~ 90°) the equation reduces to ρ = w/ 2 as expected, while in a more general case, where β = 0 and α ≫ ρ the relation reduces to ρ = (w/ 2)sinα. In principle, it is possible to determine α and β with polarisation measurements. However, in reality the duty cycle of the pulse is often too small to obtain reliable estimates (see Lorimer & Kramer 2004). Alternatively, at least for α, the relation reported by Rankin (1993) can be used: ${\mathit{w}}_{\mathrm{50}\mathit{,}\mathrm{core}}\mathrm{(}\mathrm{1}\mathrm{GHz}\mathrm{)}\mathrm{=}{\mathrm{2.45}}^{\mathrm{◦}}\mathrm{·}{\mathit{P}}^{\mathrm{-}\mathrm{0.5}\mathrm{\pm }\mathrm{0.2}}\mathit{/}\sin \mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathit{,}$(6)calculated from the observed width dependency on period for the core components of pulsars (see Sect. 5.1), which is intrinsically related to the polar cap geometry. Equation (6) is valid at 1 GHz, but can be applied at LOFAR frequencies, maintaining the same dependence, if the impact angle β ≪ ρ_core; sinα should be ignored for orthogonal rotators. Additionally, Rankin (1993), Gil et al. (1993), Kramer et al. (1994), Gould & Lyne (1998) suggested that “parallel” ρ − P relations are found if the radio emission of the pulsar can be classified and separated into emission from “inner” and “outer” cones, which seem to show different spectral properties (see Sect. 5.1 for details).

Figure 4 represents the 50% and 10% widths of the profiles as a function of the pulsar period. We show the results for LOFAR data, using only the HBA data, for which we have the largest sample. In the left panel of Fig. 4, we adopted the assumptions from Rankin (1990), later followed by Maciesiak et al. (2011), and used our interpulse pulsars (overlapping their samples of “core-single” pulsars that show interpulses) as our orthogonal rotators to calculate a minimum estimated width using a fixed dependence on the period as P^{− 1 / 2}. The interpulse pulsars in our sample are plotted in red in Fig. 4 and are shown in Fig. B.2 and labelled IP in Table B.2. The red solid and dashed lines represent the best fit of the dependency of w₅₀ and w₁₀ on P^{− 1 / 2}, which should constitute a lower limit to the distribution of pulse widths. We obtain $\begin{array}{ccc}{\mathit{w}}_{\mathrm{50}}\mathrm{(}\mathrm{150}\mathrm{MHz}\mathrm{)}& \mathrm{=}& \mathrm{(}\mathrm{3.5}\mathrm{\pm }\mathrm{0.6}{\mathrm{)}}^{\mathrm{◦}}\mathrm{·}{\mathit{P}}^{\mathrm{-}\mathrm{0.50}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.02}}\\ {\mathit{w}}_{\mathrm{10}}\mathrm{(}\mathrm{150}\mathrm{MHz}\mathrm{)}& \mathrm{=}& \mathrm{(}\mathrm{10}\mathrm{\pm }\mathrm{4}{\mathrm{)}}^{\mathrm{◦}}\mathrm{·}{\mathit{P}}^{\mathrm{-}\mathrm{0.50}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.02}}\mathit{,}\end{array}$where P is in seconds, the error is quoted at 1σ for the amplitude, and the error on the power-law index -0.50 was taken to be 0.02, following Maciesiak et al. (2011).

In the right panel of Fig. 4 we present a fit to the widths of our LOFAR sample as a function of pulse period. Because the scatter is much larger than the individual error bars, we performed a non-weighted fit. The lines represent the best fit to the data (solid line) and its 1σ dispersion (dashed lines). Here we calculated $\begin{array}{ccc}{\mathit{w}}_{\mathrm{50}}\mathrm{(}\mathrm{150}\mathrm{MHz}\mathrm{)}& \mathrm{=}& \mathrm{(}\mathrm{6.2}\mathrm{\pm }\mathrm{0.8}{\mathrm{)}}^{\mathrm{◦}}\mathrm{·}{\mathit{P}}^{\mathrm{-}\mathrm{0.1}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.4}}\\ {\mathit{w}}_{\mathrm{10}}\mathrm{(}\mathrm{150}\mathrm{MHz}\mathrm{)}& \mathrm{=}& \mathrm{(}\mathrm{16}\mathrm{\pm }\mathrm{2}{\mathrm{)}}^{\mathrm{◦}}\mathrm{·}{\mathit{P}}^{\mathrm{-}\mathrm{0.3}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.4}}\mathit{,}\end{array}$where the errors are 1σ.

The widths follow an inverse dependency with the pulsar period, consistent with previous analyses at higher frequencies (e.g. Rankin 1990; Maciesiak et al. 2011; but also Gil et al. 1993; Arzoumanian et al. 2002) and at these frequencies (Kuzmin & Losovsky 1999). In general, broadening by external effects may also be expected, even though not dominant: while we were careful to avoid evidently scattered profiles in our sample, DM smearing can also contaminate it. Finally, because our data are chosen according to detectability of the pulsars at low frequencies, a different bias in the observed sample compared with high frequencies cannot be excluded. In conclusion, our determination of the relationship between w₅₀ or w₁₀ and P is only a first step to determining the relation for the model-independent beam shape, which is to be determined when more polarisation measurements are available.

Figure 5 presents the duty cycle (w₁₀/P) of the pulsars in two bands: LOFAR HBA and L-band, for comparison, plotted against the period. The values of the duty cycle are reported in Col. 5 of Table B.2. The inverse correlation that is observed, implying that shorter period pulsars have larger beams, is also evident at LOFAR frequencies. It can be used to characterise pulsar beams and help create accurate beaming models for pulsars in the Milky Way, which would in turn constrain the Galactic population and its birth rates (Lorimer et al. 1995).

4.2. Spectral evolution of individual components

No absolute flux calibration of beam-formed LOFAR data was possible with the observing setup used for these observations. Therefore no spectral characterisation could be attempted yet. Nonetheless, we attempted a characterisation of the relative amplitudes of pulse profile components for pulsars with multiple peaks. In Table B.3 we list the pulsars for which double or multiple components can be observed and separated in at least two frequency bands. For these pulsars we selected the two most prominent peaks and calculated the evolution of the relative heights with observing frequency. We chose to select the peaks as the two most prominent maxima of the smoothed Gaussian profiles and verified by eye that we were consistently following the same peak at all frequencies. We note that as a result of subtle profile evolution (see e.g. Hassall et al. 2012), shifts of the peaks in profile longitude cannot be excluded, which we did not track here. The profile evolution is in some cases quite complex and the profiles are sometimes noisy, therefore the profiles as presented in Fig. B.1 need to be reviewed before drawing any strong conclusions based on Fig. 6. We first attempted to apply a power-law fit to calculate how the ratio between the two peak amplitudes changes with observing frequency for all the pulsars: [ P2 /P1 ] (ν) ∝ ν^γ, where P1 is the peak at the earlier phase and P2 is the peak at the later phase. This fit presented large errors, and the distribution of the spectral indices was peaked close to 0, with an average of 0 ± 2, indicating no systematic evolution despite the large scatter. A similar finding was obtained by Wang et al. (2001), who only selected conal double pulsars. They concluded that a steeper spectral index for the leading or trailing component are equally likely, arguing in favour of a same origin of the peaks in the magnetosphere, as expected if both components correspond to two sides of a conal beam. They also found a dominance of small spectral indices, with a quasi-Gaussian distribution, indicating no systematic evolutive trends. The different evolution of the peaks would then be due to geometric beaming effects. While in their case the pulsars were carefully selected so as to include only the conal doubles, in our case no such distinction was followed, so that the different relative spectral indices could also depend on a different origin of the emission regions (see also Sect. 5.3).

Given the large scatter of this result, and because our sample is heterogeneous (with double and multi-peaked pulsars), we investigated the ratios more closely. Figure 6 shows the evolution of the ratios with observing frequency for each of the pulsars used for this calculation. We ordered the pulsars into two groups, taking first the pulsars for which the profiles were already studied in previous works, and sorted by right ascension within the groups. In most cases it is apparent that the simple power law is not a good fit to the data and can be misleading if measurements are only possible at two frequencies.

Table B.3 provides a detailed summary of these measurements. We observe that the prominence of the peaks seems to shift from low to high frequencies with, in most cases, a net inversion in the dominance of P1 from LOFAR HBAs to L-band. The inversion point is also indicated in Fig. 6 by the blue horizontal line. In general, starting from LOFAR frequencies, there seems to be a trend that the peaks change from being more dissimilar amplitudes at low frequencies to becoming more similar at high frequencies. We note that Fig. 1 of Wang et al. (2001) shows that a linear trend of the peaks’ ratio with frequency fits the data well in most cases, meaning that at higher frequencies the relative amplitude of the peaks will again depart from equality. Notable changes in the observed pulse profile properties at low frequencies with respect to high frequencies have previously been observed for instance by Hassall et al. (2012), Hankins & Rankin (2010), and Izvekova et al. (1993).

An observed feature that can contribute to this behaviour was discussed by Hassall et al. (2012): they modelled the profile evolution with Gaussian components that were free to evolve longitudinally in a dynamic template. The examples presented there (two of which are also in our sample: B0329+54 and B1133+16) showed that the components change amplitudes and move relative to one another.

Hardening of the spectrum of the second peak is observed in gamma-rays in the typical case of two prominent caustic peaks (Abdo et al. 2013) and is explained with the different paths that curvature-emitting photons follow in the leading and trailing side of the profile (see e.g. Hirotani 2011), but it is not obvious that this should also follow for the radio emission.

5. Discussion

5.1. Phenomenological models for radio emission

Based on the findings discussed above, we drew some conclusions on the models that have been proposed to explain the observed properties of pulsar profiles and on some predicted effects such as radius-to-frequency mapping (RFM). These models have largely been developed based on observations performed at >200 MHz.

Rankin’s model (Rankin 1983a,b, 1986, 1990, 1993; Radhakrishnan & Rankin 1990; Mitra & Rankin 2002) proposed that the emission comes from the field lines originating at the polar caps of the pulsar, forming two concentric hollow emission cones and a central, filled, core. There is a one-to-one relation between the emission height and the observing frequency, so that at different frequencies the profile evolves, as more components come into view or disappear. Rankin based her classification on the number of peaks and polarisation of the pulse profiles. Profiles with up to five components are observed (although see Kramer 1994 for a more detailed classification). The profiles can be single (S_t or S_d based on whether the profile will become triple or double at higher or lower frequencies, respectively), double (D), triple (T or cT), tentative quadruple (cQ), and quintuple (indicated as multiple M), where c represents the presumed core origin.

Lyne & Manchester (1988) found that their data agreed with the hollow cone model, and they also observed a distinction in spectral properties between core and cone emission, or at least inner and outer emission. However, based on asymmetries of the components relative to the midpoint of the profile and the presence of so-called partial-cone profiles, they proposed that a window function defines the profile shape, and within this, the locations of emission components can be randomly distributed.

A further step in this direction was made by Karastergiou & Johnston (2007), who assumed a single hollow cone structure without core emission but instead with patches of emission from the cone rims. Emission could come from different heights at the same observing frequency, but still following RFM; the number of patches changes with the age of the pulsar: up to ~10 patches, but only at one (large) height for the young pulsars, and up to ~4–5 at ~4 different (lower) heights for the older pulsars. This also explains the narrowing of the profiles as the period increases (see Sect. 4.1), and their simulations successfully reproduce the observed number of profile components (i.e. typically N_comp ≲ 5). The central components are then simply more internal and surrounded by the external ones coming from higher up in the magnetosphere: this is why they show single peak profiles and steeper spectra. Younger pulsars have been observed to have simpler profiles, but typically with longer duty cycles than those of older ones. Karastergiou & Johnston (2007) predicted that there should be a maximum height of emission of ~1000 km. The minimum height, on the other hand, is quite varied but is close to the maximum allowed for young pulsars, which then emit only from one or two patches. Because of the width depends on the period, the opening angle of the cone would then be comparatively larger at the same height for younger than for older pulsars.

5.1.1. Radius-to-frequency mapping

In the framework of the standard models for pulsar emission, where the radio emission is predicted to come from the polar caps of the pulsar, it has been hypothesised (e.g. Komesaroff 1970 and Ruderman & Sutherland 1975) and in some cases observed (Cordes 1978) that the emission cone widens as we observe it at lower frequencies because we are probing regions further away from the stellar surface where the opening angle of the closed magnetic field lines is broader. The phenomenon is more apparent at low frequencies (<200 MHz) and is therefore ideal to study using LOFAR.

A limited observational sample has always biased the conclusions about RFM. Originally (e.g. Komesaroff 1970) it was thought that the RFM behaviour could be observed as a power-law dependence of the increase in peak separation with decreasing observing frequency, and asymptotically approaching a constant separation at high frequencies (>1 GHz). It was therefore proposed that two power laws (i.e. two different mechanisms) regulated the evolution of the pulse profile, with a break frequency at approximately 1 GHz.

Thorsett (1991) analysed this dependency using a sample of pulsars observed at various frequencies. He concluded that no break frequency seemed to be necessary to model the evolution of the pulse profile. On the other hand, a simple power law (or a quadratic dependence, indiscernible with his data) and the additional constraint of a minimal emission width (or peak separation) could fit the data at all frequencies. He obtained the following functional dependence from a phenomenological model of component separation: $\mathrm{\Delta }\mathit{\theta }\mathrm{=}\mathit{A}\mathrm{·}{\mathit{\nu }}^{\mathit{\delta }}\mathrm{+}\mathrm{\Delta }{\mathit{\theta }}_{\min }\mathit{,}$(11)where Δθ is the component separation, δ is the separation power-law index of the components, and Δθ_min the constant value at high frequency that the pulse separation tends to. The predictions for δ are quite varied depending on the theoretical model (see Table 1 in Xilouris et al. 1996).

Mitra & Rankin (2002) also did extensive work on RFM. They assumed double profiles to derive from conal emission and therefore analysed a sample of ten bright pulsars showing prominent cone components. They found that inner cones are not affected by RFM and that their component separation does not vary with observing frequency, while outer cone components show RFM and increase their separation with decreasing observing frequency.

Hassall et al. (2012) have discussed that RFM does not seem to be at play for some pulsars observed with LOFAR. Here, with a more statistically significant sample, we can discuss the matter in more detail. Following Thorsett (1991) and in particular Mitra & Rankin (2002), we investigated double-peaked pulse profiles and their component separation (see peak phases in Cols. 5 and 6 of Table B.3).

Mitra & Rankin (2002) divided a group of double profile pulsars into three groups The pulsars from groups A and B are associated with outer cone emission (the difference between the two being a fit with or without the constraint ρ₀ = ρ_pc, where ρ₀ is the constant equivalent to Δθ_min from Eq. 11, relative to the beam radius, and ρ_pc is the beam radius at the polar cap edge) while the pulsars from group C are associated with inner cone emission. Of the ten pulsars of Mitra & Rankin (2002), the pulsars from our sample that fall inside each group are

Group A: B0301+19, B0525+21, and B1237+25

Group B: B0329+54 and B1133+16

Group C: B0834+06, B1604−00, and B1919+21.

Mitra & Rankin (2002) reported that pulsars from groups A and B show RFM, while the pulsars from group C show almost no evolution at all. Although the profile in our sample evolves rapidly at low frequency, it seems that a similar behaviour can be observed (see single cases in Fig. B.1).

Figures 3 and 7 show a similar calculation using w₁₀. We plotted the evolution of w₁₀ as a function of observing frequency for the single-peaked pulsars and the histogram of the spectral indices of this evolution. We excluded the LBA and HBA profiles that showed significant scattering. The errors were calculated from the Gaussian fit, taking into account both the noise contribution and any unaccounted scattering of the profile. Although the values in Fig. 3 seem to follow the power-law, the profile width in some cases effectively behaves in a non-linear way, as can be cross checked in Fig. B.1 for the single cases.

The weighted mean spectral index from Fig. 7 is δ = −0.1(2). Our results are compatible at 1σ with the predictions made by Barnard & Arons (1986): the component separation does not vary (no RFM, δ ~ 0.0), but the distribution in Fig. 7 peaks at negative spectral indices, which is evidence for a weak widening of the profile at low frequencies. Following the predictions of Gil & Krawczyk (1996), the dependence of w₁₀ with observing frequency based on their calculations should be δ = −0.21 for RFM and conal beams. This is in support of their model, while the fact that we see a broad distribution and a flatter median index might be explained by a subdivision of pulsar behaviours according to the Rankin groups. For a future complete analysis, geometry should be taken into account to perform a study on the beam radii (ρ) rather than the pulse widths.

An alternative but complementary explanation to the observed widening of pulsar profiles with decreasing observing frequency can be found in the theory of birefringence of two different propagation modes of a magneto-active plasma (McKinnon 1997). These two modes of propagation follow a different path along the open field lines. The nature of birefringence is such that the two polarisations are spatially closer together at higher frequencies, and depolarisation will occur where they overlap. Beskin et al. (1988) predicted that the two modes of propagation would result in two different indices: δ = −0.14 and δ = −0.29 for the ordinary mode and δ = −0.5 for the extraordinary mode. While our observations are compatible with either scenario or a combination of them, polarisation studies will help discern between the two interpretations of this phenomenon (see Noutsos et al. 2015).

5.1.2. Profile complexity

Karastergiou & Johnston (2007) searched for a relation between the number of peaks in the profile of a pulsar and some observed or derived parameters, such as its period, age, and rotational energy loss. As a general trend, they found that faster, younger, more energetic pulsars would typically show less complex profiles, which prompted them to assume that the regions of emission for these pulsars arise at higher altitudes in the magnetosphere and are, therefore, less numerous (see Sect. 5.1). Notably, an abrupt change in this respect can be observed at P< 150 ms, τ< 10⁵ yr and Ė> 10³⁵ erg s^-1.

We calculated the same relations using our sample of 100 pulsars and made the comparisons using LOFAR HBA band and the L-band data. Figure 8 shows the relation between the period of the pulsar and its spin-down age τ for the two frequency bands. Each circle represents a pulsar, and its colour and diameter represent a different number of peaks of its profile. The circles are larger for growing number of peaks with frequency, and the number of peaks is ordered by colour, in the order white, red, green, blue, and purple. In the histograms we summed the pulsars separated by number of peaks in their profiles to search for trends as a function of either period or τ. No trends are evident in any of the histograms. Our sample does not cover the region of young energetic pulsars in a statistically significant way so that while the few cases might confirm the predictions of Karastergiou & Johnston (2007), nothing in favour or disfavour can be stated in this respect.

Comparing the histograms from the two plots it appears, in general, that the number of peaks changes from one band to the other. For more clarity, we plotted the trend of the number of peaks from HBAs to L-band in Fig. 9. Here the triangles pointing down (red) represent the pulsars that show fewer peaks at HBA frequencies than at L-band: a total of 30 pulsars; and the ones pointing up (blue) are those for which the number of peaks is higher at HBA frequencies than at L-band: 16 pulsars. The green triangles pointing right represent the pulsars whose number of peaks does not change between the two frequencies: 30 pulsars. There is a slight indication that faster-spinning pulsars have a higher number of peaks at L-band than at HBA frequencies while the slower spinners have more peaks at HBA frequencies than at L-band. The younger pulsars in general have more peaks at L-band than HBAs; the older have fewer peaks at L-band than at HBAs. In both cases the null hypothesis probability that the two distributions are the same based on the Kolmogorov-Smirnoff statistics cannot be rejected at the 20% confidence level. This, as pointed out by Lyne & Manchester (1988), might be a selection effect or a resolution effect (as fast pulsars tend to have wider pulses that can be resolved more easily), but given our statistics, we can also tentatively assume that our sample includes pulsars with different behaviours.

Additionally, a number of pulsars show an increase in the number of components at HBA compared to L-band, which may also be explained at HBA frequencies by the general expectation that a higher portion of the beam can come into view at lower frequencies, according to RFM. On the other hand, the fastest and youngest pulsars show an opposite trend, similar to what is also observed in millisecond pulsars (see Kondratiev et al. 2016), which are even faster and therefore have a wider duty-cycle (see Fig. 5). A similar finding was also reported by Hankins & Rickett (1986), who analysed the frequency dependence of pulsar profiles for 12 pulsars in the 135–2380 MHz frequency range. They argued that the occurrence of single profiles at low frequencies that become multiple at high frequencies can be explained in the framework of the “core and cones” models (in the formulation of Rankin 1993). In particular, this is expected to occur if we can only observe the core emission at low frequencies, which is observed to typically have a steeper spectrum, while the outriding conal components only emerge at higher observing frequency (see Kramer et al. 1994 for arguments why this is caused by geometrical reasons and applies to inner and outer components regardless of their “nature” as core or cones). A word of caution is needed here, related to how the number of peaks was determined: it is possible that we achieved a good fit to the profile using a smaller number of Gaussians in HBAs relative to L-band because the quality of the profile is lower and so fewer components need be fit (for details on the method see Sect. 3, and the single cases can be studied by comparing Fig. B.1 and Table B.2).

5.2. External effects on profile evolution

The interstellar medium affects the pulse signal while it travels towards the observer, and it strongly depends on observing frequency, with observations at low frequencies being more strongly affected by scattering and dispersion delay (see also Zakharenko et al. 2013). We here did not correct the profiles presented for scattering effects that can smear out the signal especially at LBA frequencies, but we considered the effect of intra-channel dispersive smearing.

The DM represents the integrated column density of free electrons between the pulsar and the observer. It produces a time delay in the signal, between the observing frequency and infinite frequency, that can be approximated as $\mathrm{\Delta }{\mathit{t}}_{\mathrm{DM}}\mathrm{=}\frac{\mathrm{\left[}\frac{\mathrm{DM}}{\mathrm{c}{\mathrm{m}}^{-3}\mathrm{pc}}\mathrm{\right]}}{\mathrm{2.41}\mathrm{\times }{\mathrm{10}}^{-4}\mathrm{[}\frac{\mathit{\nu }}{\mathrm{MHz}}{\mathrm{]}}^{\mathrm{2}}}\mathrm{s}\mathit{.}$(12)This approximation is valid if the plasma is tenuous and thus collisionless and if the observing frequency is much greater than the plasma frequency and the electron gyrofrequency. Signals at different frequencies will be delayed, with the lowest frequencies being delayed the most. These delays can change on timescales of some years, up to 10^-3 cm^-3 pc (see e.g. Keith et al. 2013).

When aligning the profiles absolutely as we did, connecting the reference point of the profile to the reference epoch of the ephemeris, the DM delays had to be taken into account and all reference times converted to the corresponding times at infinite frequency to correct for dispersive delay. Hankins & Rickett (1986) described a way to measure DM variations based on the alignment of pulsar profiles at different frequencies. Hankins et al. (1991) and Hankins & Rankin (2010) followed this method using increasingly higher resolution profiles. They absolutely aligned the profiles spanning, where possible, all the octaves of radio frequency. Leaving the DM as the only free parameter, they identified a reference point in the profile and adjusted the DM value to compensate for the remaining misalignment of the profiles. The value of the best DM was obtained this way, its accuracy strongly depending on the lower frequency that can be used and on the precision with which the time difference of the misalignment can be measured.

All the pulsars from our sample were aligned by refolding their profiles at all frequencies using the same ephemeris. The DM that was adopted for the alignment was the one obtained from LOFAR HBA data: a first folding was performed with prepfold using the new ephemeris created from the Lovell data (see Sect. 3), but allowing for a search over DM values, and then the best DM obtained from this search was included in the ephemeris and the profile was dedeispersed once more, without any search option. The same was done for the LBA and the P- and L-band observations. Figure B1 shows that the multi-frequency profiles are aligned in most cases. There were some cases, discussed in Sect. 3, where the alignment is visibly incorrect or where we had to apply an extra adjustment to DM to compensate for a visible offset between the profiles (see Table 1).

One possible cause for the misalignment is the evolution of the profile across the frequencies, so that it is not possible to easily identify a fiducial point in the profile. In addition, pulsars strongly affected by scattering will not only have an exponential scattering tail, but also an absolute delay. Our profile alignment seemed to also be affected by the significantly different observing epochs (e.g. the WSRT observations were performed more than ten years before LOFAR observations). On one hand, this meant that we had to adopt a timing solution spanning a long period of time, where timing noise and other effects can become substantial. On the other hand, with observations so far apart, we might also be probing gradients in electron density.

In our case it is not yet possible to perform a systematic study of these variations, as was done by Keith et al. (2013). Nonetheless, LOFAR data can provide a new wealth of DM measurements to be compared with previous observations to map the evolution of the interstellar dispersion with time. A first conclusion that can be drawn from this, simply by comparing the DM values in Table B.1, is that there is no significant indication of DMs being systematically higher at low frequencies (at least at our measurement precision). Some authors (e.g. Shitov 1983) have postulated “superdispersion” due to the sweepback of field lines in the pulsar magnetosphere, which would be responsible for lower dispersion delays at low frequencies and could create an observed profile misalignment over a wide observing band. We find that the ratio DM_HBA/ DM_eph ranges from 0.97 to 1.06 (see Table B.1), thus differing by a significant percentage in some cases, but in both directions, thus not favouring superdispersion. This agrees with previous findings using LOFAR data (Hassall et al. 2012) and previous measurements (Hankins et al. 1991).

5.3. Some examples of unexpected profile evolution

It is not in the scope of this initial paper to enter into much detail about the profile evolution of specific pulsars. These will be the subject of future dedicated work. There are, however, several interesting cases of pulse evolution that are worth pointing out at this time.

Figures 3 and 7 showed that there are some cases (e.g. B1541+09, B1821+05, and B1822–09, B2224+65) where the width of the profile is observed to increase with increasing observing frequency. If we compare these results with the single profiles in Fig. B.1, we notice that in these cases this is caused by new peaks appearing in the profile at higher frequencies (e.g. the well-known “precursor” in the case of B1822–09). This is not common in the standard core and cones geometries, even though it is predicted that new components can come into view if our viewing angle changes with increasing observing frequency, thus allowing us to see deeper into the beam. While this would explain some of the cases, in some others the profile evolution can hardly be ascribed to a symmetric core and conal structure (see e.g. B0355+54, B0450+55, B1831–04, and B1857–26).

These narrowing profiles at low frequencies might be interpreted as evidence for fan beam models (Michel 1987). In fan beam models by Dyks et al. (2010), Dyks & Rudak (2012, 2015), for instance, the emission comes from elongated broad-band streams that follow the magnetic field lines. The model, based on the cut angle at which the line of sight crosses the beam, can explain the lack of RFM, for example, if the stream is very narrow (also the case for millisecond pulsars), and it can also explain the “inverse” RFM if there is spectral non-uniformity along the azimuthal direction of the beam through which our line of sight cuts. The fan-beam formulation proposed by Wang et al. (2014) can even explain “regular” RFM by assuming that a fan beam composed of a (small) number of sub-beams will produce a so-called “limb-darkening pattern” that is caused by the decrease in intensity with beam radius of the emission at higher altitudes because the emission is farther from the magnetic pole. Their model, based on observations and simulations, predicts that the non-circularly bound beam (different in this from the beam predicted by the narrow-band models) can depart from the relation w ∝ P^{− 1 / 2} (where w is the measured width of the profile, see Sect. 4.1).

Chen & Wang (2014), who recently analysed the pulse width evolution with frequency of 150 pulsars from the EPN database, reached a similar conclusion: the emission must be broad-band, and the observed behaviour of width at different frequencies is caused by spectral changes along the flux tube. They found that if the spectral index variation along pulse phase is symmetric, there can be either canonical RFM or anti-RFM, while in cases where there are substantial deviations from the symmetric case, then there can be the non-monotonic trends of w₁₀ with ν, which we also observed. This is supported theoretically by the particle-in-cell simulations of pair production in the vacuum gaps (Timokhin 2010), which predict that the secondary plasma does not necessarily have a monotonic momentum spectrum.

These results, and in particular the fact that a wide stream is expected to produce spectral variations longitudinally in the profiles, could also explain the observed behaviour of the peak ratios shown in Fig. 6. Alternatively, the peak ratio changing with observing frequency and, in particular, its changing sign, might be related to the frequency dependance of the two modes of polarisation (ordinary, “O”, and extraordinary, “X”) (e.g., Smits et al. 2006) and to that they might be differently dominant in different peaks. While it is not possible to give a comprehensive analysis of the phenomenon here, our studies on the polarised emission from pulsars with LOFAR (see Noutsos et al. 2015) address the questions related to the orthogonally polarised modes and the related jumps in the polarisation angle.

While in this section we have discussed some unexpected profile evolution, it remains the case that most pulsar profiles (cf. Sects. 4.1 and 5.1.1) are well described by RFM.

6. Summary and future work

We have presented the profiles of the first 100 pulsars observed by LOFAR in the frequency range 119–167 MHz. Twenty-six of them were also detected with 57 min integrations using LOFAR in the interval 15–63 MHz. All the LOFAR profiles presented in this work will be made available through the EPN database¹⁰. LOFAR observations were compared with archival WSRT or Lovell observations in P- and L-band, after first folding and aligning all profiles using an ephemeris spanning the full range of the observations. The rotational and derived parameters are presented in Table B.1. Two values of DM were presented as well: one obtained by the best timing fit and one from the best fit of LOFAR data. For each pulsar we aligned the profiles at different frequencies in absolute phase, using the latter DM value. The 100 profiles are presented in Fig. B.1.

Each profile of every pulsar was described using a multi-Gaussian fit following the approach of Kramer et al. (1994), so that in general more components were needed to fit the profiles than evident at first glance or traditionally considered (e.g. Rankin 1983b). The results of the Gaussian fit (the measurement of the widths at half and at 10% of the maximum of the profile and the spectral index of their evolution with observing frequency) are reported in Table B.2. Using the components’ widths, we calculated the ratio of the peaks for pulsars with double or multiple peaks (the two most prominent, in the case of multiple peaks pulsars). We concluded that the ratio of the main peak to the second peak does not follow a unique trend, although we note that in most cases the dominant peak alternates with changing observing frequency. Using w₁₀ , we followed the evolution of the width of the full profile with observing frequency. We concluded that while our average spectral index is compatible with no evolution of the pulse width, the distribution of the values is quite large and compatible with the presence of different behaviours for different pulsars, for example, based on inner or outer cone emission, as discussed by Mitra & Rankin (2002), or on different propagation modes in the magnetosphere (e.g. Beskin et al. 1988; Beskin & Philippov 2012).

Future work will be needed, and is in progress, to add more elements to complete this puzzle. Parallel to this work, a similar one is being conducted on the evolution of the profiles of millisecond pulsars at low frequencies (Kondratiev et al. 2016), and the spectral behaviour of the slow and recycled pulsars has been analysed (Hassall et al., in prep.). Complementary to our work is the study of the polarisation properties of pulsars (Noutsos et al. 2015): the study of the polarisation properties can give constraints on the geometry of the pulsar emission and therefore on its height and on the intrinsic opening angle of the beam of the emission. Additionally, polarisation can help distinguish between orthogonal polarisation modes and therefore determine whether the observed widening of the profiles is caused by birefringence.

To better constrain the width of the pulse, it is also important to be able to deconvolve the scattering tail, which in some cases becomes dominant at low frequencies, from the intrinsic width of the profile. Studies on the characterisation of scattering at low frequencies and modelling of the scattering tail are being conducted (Archibald et al. 2014, Zagkouris et al. in prep.). At the same time, the effects of the ISM on LOFAR profiles are being used to create an “ISM weather” database, where the DM variations, independently measured with LOFAR (Verbiest et al., in prep.), can be used for high-precision timing measurements from pulsar timing arrays (see Keith et al. 2013).

Finally, the observations presented in this work were from commissioning LOFAR data; as mentioned, LOFAR is currently performing the LOTAAS all-sky pulsar survey. When completed, it will provide an extraordinary database of 1.5 h coverage of the whole northern sky with 0.49 ms time resolution and down to a flux S_min ~ 6 mJy at 135 MHz. At present we are able to use the full LOFAR core instead of only the “Superterp” and can cover a full 80 MHz bandwidth contiguously. Moreover, further improvements will include coherent dedispersion of the data and will enable single pulse studies to probe the “instantaneous” magnetosphere.

---

Acknowledgments

M.P. wishes to thank A. Archibald, V. Beskin and J. Dyks for useful discussion. We are grateful to an anonymous referee whose comments notably improved the quality of our work. This work was made possible by an NWO Dynamisering grant to ASTRON, with additional contributions from European Commission grant FP7-PEOPLE-2007-4-3-IRG-224838 to JVL. 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. M.P. acknowledges financial support from the RAS, Autonomous Region of Sardinia. J.W.T.H., V.I.K. and J.V.L. acknowledge support from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement nrs. 337062 (JWTH, VIK) and 617199 (JVL). SO is supported by the Alexander von Humboldt Foundation. C. Ferrari acknowledges financial support by the “Agence Nationale de la Recherche” through grant ANR-09-JCJC-0001-01.

References

Appendix A: Determination of pulse widths

As discussed in the Sects. 3 and 4.1, the widths of the full profile that are to be used in the subsequent calculations were determined from the fit of Gaussian components to the profile shape. However, a number of alternative methods were also tested and, in one case, used to cross check the widths obtained from the Gaussian fits. Here we give a brief description of each of them.

We calculated the effective width, w_eff, as the integrated pulsed flux divided by the peak flux; this is represented by the cyan-shaded areas in Fig. A.1 centred at the main peak. This width metric typically does not represent all the profile components well, as is observed in the P- and L-band profiles. The method underestimates the width of the profile for multi-peaked profiles, which is the majority of profiles in this work, mainly characterising the width of the main peak in these cases.

The total-power width w_pow is calculated as the phase interval that includes 90% of the total pulse energy in the cumulative flux distribution of the normalised profile (see lower plots in Fig. A.1). We selected the phase included in the interval between 5% of the total flux and 95% of the total flux (indicated by the intersection of the summed profile with the horizontal lines in panel b of Fig. A.1). This distribution should start at 0

and increase monotomically to 1 if no noise is present (see the noise-free distribution on the right side of Fig. A.1). When noise is present (left-hand side of Fig. A.1), the sum does not grow smoothly, and sometimes negative terms in the off-pulse region add up quite substantially. In the real profiles (left-hand side of Fig. A.1), the horizontal lines should all overlap at 0.05 and 0.95, like they do in the case of the noise-free profiles (right-hand side of Fig. A.1), but they differ in some cases as they take into account that the highest and lowest value of a noisy cumulative distribution can be negative or greater than 1. As can be seen from the green (LBA) and magenta (P-band) curves, when the profile is noisier, this measurement is less stable because the cumulative distribution oscillates more, and the width can be consistently overestimated (in these cases the left-hand dashed line precedes the phase range shown in the Fig.). This method is represented by the dotted line in the upper plot of Fig. A.1, which, in the case of the LBA and P-band profiles, starts at phase 0.

The full-width at 10% of the maximum, w_op, represents the width of the full profile (including noise) at the 10% level of the outer components, including the full on-pulse region: it is coincident with the solid vertical lines in Fig. A.1. It was calculated as the on-pulse region with a flux above 10% of the main peak in the baseline-subtracted profile.

Appendix B: Multi-frequency profiles and tables of the derived profile properties.

All Tables

All Figures