Mode switching and oscillations in PSR B1828-11

TL;DR

This study reveals that PSR B1828-11 exhibits mode switching between two stable states, with transition rates linked to spin-down variations, suggesting magnetospheric changes as the cause rather than free precession.

Contribution

The paper provides detailed multi-hour observations demonstrating mode switching and links it to spin-down changes, challenging previous free precession models.

Findings

Mode switching occurs between two stable profile states.

Transition rates correlate with spin-down variations.

Free precession is not required to explain the observed phenomena.

Abstract

The young pulsar PSR B1828-11 has long been known to show correlated shape and spin-down changes with timescales of roughly 500 and 250 days, perhaps associated with large-scale magnetospheric switching. Here we present multi-hour observations with the Parkes and Green Bank Telescopes at multiple phases across the roughly 500-day cycle and show that the pulsar undergoes mode-changing between two stable, extreme profile states. The fraction of time spent in each profile state naturally accounts for the observed overall "shape parameter" (defined to be 0 for wide profiles and 1 for narrow ones); this and the variable rate of the mode transitions are directly related to the spin-down changes. We observe that the mode transition rate could plausibly function as an additional parameter governing the chaotic behaviour in this object which was proposed earlier by Seymour and Lorimer. Free precession is not needed to account for the variations.

Tables

Table 1. Table 1: Observations and mode-switching properties. Uncertainties in the last quoted digit(s) are given in parentheses. The entries in the “Valid Data” column are smaller than those in the “Length” column for those epochs in which the data are split into several short files. Because dspsr labels the integrations with integer multiples of 10 seconds, the data length values are generally very slight overestimates.

MJD	Telescope	$\dot{P}$ cycle	Length	Valid	Average	Number of	Transition	Peak ratio
		phase	(s)	data (s)	state	transitions	rate (s-1)	(W/N)
51579	Parkes	0.330(5)	1510	1510	0.901	14	0.00927	0.74
51588	Parkes	0.348(5)	4840	4520	0.781	83	0.01836	0.64
51633	Parkes	0.440(5)	11920	10610	0.445	45	0.00424	0.48
51742	Parkes	0.662(5)	3320	2980	1.000	0	0.00000	0.0
51757	Parkes	0.693(5)	10180	9050	1.000	0	0.00000	0.0
51842	Parkes	0.867(5)	14000	12000	0.021	2	0.00017	0.46
52041	Parkes	0.275(5)	6620	6010	1.000	0	0.00000	0.0
52132	Parkes	0.462(5)	13560	11140	0.458	36	0.00323	0.45
52189	Parkes	0.580(5)	15610	13180	0.351	10	0.00076	0.46
52261	Parkes	0.728(5)	14400	14400	1.000	0	0.00000	0.0
52318	Parkes	0.846(5)	13070	13070	0.000	0	0.00000	—
52394	Parkes	0.003(5)	14060	14060	0.000	0	0.00000	—
52408	Parkes	0.033(5)	13950	13950	0.000	0	0.00000	—
52466	Parkes	0.152(5)	14400	14400	0.287	1	0.00007	0.46
52532	Parkes	0.289(6)	10610	10610	1.000	0	0.00000	0.0
52625	Parkes	0.483(6)	1350	1350	0.533	4	0.00296	0.48
53989	GBT	0.358(10)	15650	15650	0.473	152	0.00971	0.54
54047	GBT	0.482(11)	15900	15900	0.379	29	0.00182	0.46
54127	GBT	0.653(11)	15270	15270	1.000	0	0.00000	0.0
54202	GBT	0.814(12)	17050	17050	0.119	12	0.00070	0.45

Table 2. Table 2: Flux densities for the GBT data. Uncertainties are estimated to be about 10%.

MJD	Wide	Narrow	Total	Ratio
	(mJy)	(mJy)	(mJy)	(W/N)
53989	1.3	1.5	1.4	0.85
54047	1.0	1.3	1.13	0.74
54127	—	1.5	1.5	—
54202	1.07	1.5	1.12	0.73

Table 3. Table 4: Dominant significant periodicities and corresponding timescales and numbers of pulse periods found in epochs with multiple transitions.

MJD	Frequency	Timescale	Periods
	(Hz)	(s)
51633	0.0002435	4108	10140
	0.0003610	2770	6340
	0.0004618	2165	5345
52132	0.0004649	2151	5311
52189	0.0000833	12005	29641
	0.0002372	4216	10410
53989	0.0005179	1931	4768
54047	0.0002958	3381	8347
	0.0004091	2444	6036
	0.0005223	1915	4727
54202	0.0003228	3098	7649

Table 4. Table 5: Fit to period derivative data. Uncertainties in the last quoted digit are given in parentheses.

Parameter	Value
A	59.993(2) $\times {10}^{-15}$ s/s
B	$-1.6$(2) $\times {10}^{-15}$ s-1
C	$-0.118$(3) $\times {10}^{-15}$ s/s
D	$-0.102$(3) $\times {10}^{-15}$ s/s
$P_s$	482.9(9) days
${\dot{P}}_s$	$-0.0048$(5) days/day
$t_1$	MJD 51907(2)
${\varphi }_2$	0.29(5) rad

Equations

$$\dot{P}=A+B(t-52930)+C\cos \left(\frac{2\pi (t-t_1)}{P_s+\dot{P_s}\times (t-52930)}\right)$$

$$+D\cos \left(\frac{2\pi (t-t_1)}{(P_s+\dot{P_s}\times (t-52930))/2}+{\varphi }_2\right),$$

Full text

††pubyear: 2017††pagerange: Mode switching and oscillations in PSR B1828$-$11–Mode switching and oscillations in PSR B1828$-$11

1 Introduction

A prime characteristic of radio pulsars, and the property that makes them valuable tools for the precision study of physics ranging from the neutron-star interior to cosmological structure formation, is the stability of their rotation periods. Pulsar timing studies enumerate every rotation of a neutron star across years or decades with no ambiguity; this leads directly to high-precision measurements of spin periods, astrometric parameters and binary parameters, if applicable. Radio pulsars also generally have extremely stable pulse profiles at any given observing frequency. This property is exploited to provide well-measured pulse arrival times by cross-correlation with a low-noise “standard profile”.

There are some notable exceptions to the pulse stability rule, particularly when individual pulses are observed. Some pulsars “null," disappearing for anywhere from a few pulses to weeks at a time (e.g. Backer, 1970b; Kramer et al., 2006). Another common behaviour is that of mode-changing, in which a pulsar displays two or more “typical" profiles and switches between them with varying timescales (Backer, 1970c). There is an indication of a relationship between these phenomena (e.g., van Leeuwen et al., 2002) and also between the drifting subpulses that sum together to produce average profiles in a number of bright pulsars.

The young pulsar PSR B1828$-$11 (PSR J1830$-$1059) was discovered in a 20-cm survey of the Galactic plane conducted with the 76-m Lovell Telescope at Jodrell Bank Observatory in the late 1980s (Clifton et al., 1992). Its pulse period is 405 ms and its dispersion measure is close to $160$  pc cm*-3*. PSR B1828$-$11 has a characteristic age of 110 kyr and an inferred surface magnetic field strength of $5\times 10^{12}$ G, making it at first glance an unremarkable object among the pulsar population. However, routine and eventually more intensive follow-up with the Lovell Telescope showed that it displays roughly periodic variations in its spin-down rate (period derivative), with periodicities of roughly 250 and 500 days. These were correlated with the observed pulse shape, which ranged between “wide" and “narrow" extremes on the same timescales (Stairs et al., 2000). The initially favoured interpretation was that free precession was causing the spin axis to wobble, providing variable torque on the neutron star and allowing the observer to see different parts of the emission region over time. This idea was not without problems, notably the expectation that pinned vortices in the superfluid internal to the neutron star would damp out such precession on timescales of several hundred precession cycles (e.g., Shaham, 1977; Sedrakian et al., 1999). Later authors (e.g. Jones & Anderson, 2001; Link & Epstein, 2001; Link & Cutler, 2002; Link, 2006) argued that the existence of precession should instead be viewed as setting a limit on the amount of pinned superfluid in the star, with possibly large implications for neutron-star structure. A complicating factor is that this pulsar underwent a glitch on 2009 July 29 (MJD 55041; Espinoza et al., 2011, and http://www.jb.man.ac.uk/~pulsar/glitches/gTable.html), which may force a requirement for pinned vortices after all (Jones et al., 2017).

Meanwhile, other pulsars were beginning to show signs of correlation between profile shape and spin-down properties. The most notable of these is PSR B1931+24, which nulls for weeks on end, and moreover has a spin-down rate 50% larger when the pulsar is on (Kramer et al., 2006). Further examples of this intermittent behaviour are now known in four other pulsars (Camilo et al., 2012; Lorimer et al., 2012; Lyne et al., 2017). In the case of PSR B1931+24, it seems clear that global changes in the magnetosphere affect both the torque on the neutron star and the ability of the pulsar to produce radio emission. Using the Lovell Telescope pulsar data archive, Lyne et al. (2010) argued that a similar type of magnetospheric switching must be taking place for six other pulsars including PSR B1828$-$11. This paper also presented evidence, from the average profiles of short observations, that PSR B1828$-$11 likely spends most of its time in the extreme wide or narrow states, with successive days sometimes producing different profiles.

Another model put forth to explain the phenomena seen in PSR B1828$-$11 involves non-radial oscillations (Rosen & Clemens, 2008; Rosen et al., 2011) in which pulsations of the neutron star with high spherical degree can modulate the emission and produce different pulse modes and/or drifting subpulses while simultaneously changing the torque.

More recent work has compared the Bayes factors for a precession model and a magnetospheric switching model for PSR B1828$-$11, not taking into account any short-term changes and finding a weak preference for the precession model (Ashton et al., 2016). Ashton et al. (2017) further show that the precession period is decreasing, which would imply an increase in the stellar deformation. Jones (2012) argues that state switching and precession are fully compatible, with the latter driving changes in the time spent in each magnetospheric state. Kerr et al. (2016) similarly advocate that PSR B1828$-$11 is undergoing precession which in turn is regulating the magnetospheric switching.

Here we report on multi-epoch long observations of PSR B1828$-$11 with the Parkes and Robert C. Byrd Green Bank Telescopes which present evidence for mode-changing which differs from epoch to epoch, relate the mode-changing to the timing variability and show that the chaos model may be a good descriptor of the data. Preliminary analyses of some of these data were presented in Stairs et al. (2003) and Lyne (2013) and a brief overview of these results can be found in Stairs et al. (2018). In §2, we describe the observations and analysis. In §3, we present our results which are then discussed in detail in §4. Finally, in §5, we draw conclusions and suggest directions for future work.

2 Observations and Data Reduction

At the 64-m Parkes telescope, data were taken at 1390 MHz between 2000 Feb. 5 and 2002 Dec. 17 with a $2\times 512\times 0.5\,$MHz filterbank which provided 1-bit samples every 0.25 ms. The two polarizations were summed before digitization and the data are not flux-calibrated. Most of these observations are about 4 hours in length, with the first ones broken up into several short files, while the data spans on a few days are shorter. At the 100-m Green Bank Telescope (GBT), data were collected between 2006 Sept. 11 and 2007 April 12 with the Berkeley-Caltech Pulsar Machine (BCPM; Backer et al., 1997) with 96 1-MHz frequency channels centred at 1400 MHz. Polarizations were summed in hardware and the data stream was sampled every 72 $\mu s$. All these observations were about 4 hours in length. Accompanying short calibration files were taken using a 25 Hz switched diode with temperature 1.6 K. Uncalibrated data were used to determine the switching states and derived information, but flux-calibrated summed profiles and the resulting flux densities are presented later in this work.

The observation dates were chosen to sample the roughly 500-day period-derivative cycle. In Fig. 1 we plot the observing epochs over the timing and average-pulse-width data presented in Lyne et al. (2010). Observations and basic properties are summarized in Table 1, where we also include the phase of each observation relative to a fit to the $\sim 500$-day $\dot{P}$ cycle (see Section 3.2) as computed for the integer MJDs.

For each epoch, the filterbank search-mode data were folded using dspsr (van Straten & Bailes, 2011) with 10 s sub-integrations and a dispersion measure of $159.7$  pc cm*-3* (Stairs et al., 2000). This sub-integration length allowed the straightforward identification by eye of “wide" and “narrow" states within each observation using PSRCHIVE’s pav and paz programs (van Straten et al., 2012). An example plot of a short Parkes data file is shown in Fig. 2; the times at which the mode changes occur are just as clear in the other data sets.

3 Results

3.1 Mode changes

In Fig. 3 we show the switches between wide (0) and narrow (1) modes at each epoch. Breaks in the plots occur when an observation is broken up into multiple data files. It is clear that the fraction of time spent in each mode varies from epoch to epoch, along with the frequency of transitions. Table 1 includes the average mode value along with the transition rate.

At each epoch, the data files were subdivided into narrow and wide subfiles, and then summed to produce average profiles. These are shown in Fig. 4. Although the Parkes data are not flux-calibrated (the levels are roughly consistent from epoch to epoch), it appears that the overall pulse heights are reasonably stable. Table 2 presents the flux densities for the cumulative wide and narrow profiles for each GBT epoch. The two days with long-duration sequences (see Fig. 3) have similar flux density ratios between the wide and narrow pulses. The small overall difference in flux density may be due to refractive scintillation. The epoch with the most rapid transitions, MJD 53989 (2006-Sept.-11), has a larger relative flux density in the wide profile; this is consistent with the higher peak of wide profile in Fig. 4. Note that a similarly higher-peaked wide profile occurs for the Parkes epoch with the most rapid transitions, namely MJD 51588 (2000-Feb. 14). Whether or not this represents a true change in the wide-mode profile is not immediately clear, given the presence of significant radio-frequency interference on MJD 53989 and the inevitable inclusion of some narrow pulses within the wide-designated sub-integrations due to the 10-second sub-integration length. This latter systematic would have the effect of increasing the height of the wide profile while simultaneously reducing the height of the narrow one. The apparent small shape changes in the wide component may be due to a combination of noise and, again, the admixture of the narrow pulses contaminating some sub-integrations. The small trailing “shoulders” on the narrow components on these two days also suggest such cross-contamination. In general, it is not possible to distinguish single pulses in these data (but see § 4.5), and in any case we are most interested in the pulse envelope. The narrow and wide pulse profiles look very similar in shape at all epochs.

To quantify this similarity, we have determined how well the per-epoch average profiles are fit by linear combinations of the extreme wide and narrow states, choosing the pairs of profiles on MJD 52466 for the Parkes epochs and on MJD 54204 for the GBT epochs as the standard profiles, as these dates had good representation of each state, with only a few, completely unambiguous changes. The resulting profile residuals are shown in Fig. 5. Their small relative amplitudes demonstrate that the two extreme profiles are sufficient to model the pulse. This leads us to conclude that there are most likely only two profile envelopes visible to us no matter what the $\dot{P}$ phase. This argues against a precession interpretation in which the impact angle to the magnetic pole is smoothly varying.

Histograms of the wide and narrow segment durations derived from the fully-observed wide and narrow segments (i.e., those not at the boundaries of the observation) seen for each epoch with significant switching are shown in Fig. 6. With these we plot the probability density functions (with arbitrary scaling) for two-state Markov processes with means derived from the histograms, as described in Cordes (2013). In three cases (the wide states on MJDs 51579, 51588 and 53989 – all at roughly the same phase in the $\dot{P}$ cycle), these curves appear to be reasonable matches to the distributions, but on most days that is not the case. The histograms do not appear to be well fit in general by exponential distributions, as would be expected if the mode transitions followed Poisson statistics. There is often a pattern in which the wide sequences are more prone to very short durations, but this is not universally true.

We carried out 2-sided Kolmogorov-Smirnov tests using the scipy routine stats.ks_2samp, to check whether the complete wide and narrow segments were drawn from different distributions, and whether the narrow and wide segment distributions differ from epoch to epoch. We summarize the results in Table 3. Overall there is evidence for some epochs having similar state durations and some having very different ones, as might be expected from the data in Fig. 3. When considering these results, it should be borne in mind that MJD 51579 has very few full segments (7 wide segments and 6 narrow).

In an effort to explore the long-period spectral content within the observations, we computed discrete Fourier Transforms of the state sequences for each epoch with rapid changes along with Lomb-Scargle periodograms using the astropy routine LombScargle (Astropy Collaboration: Robitaille et al., 2013). These are shown in Fig. 7. While the multiple short observations employed in the earliest Parkes data sets result in some unfortunate windowing (cyan curves), it is nevertheless clear that there are significant periodicities of the order of 30 minutes to an hour in several of the data sets. In Table 4 we provide the frequencies and corresponding numbers of pulse periods for the most dominant well-defined peaks with false-alarm probabilities of less than 0.0001 in the Lomb-Scargle analyses, leaving out those which coincide with the relevant window functions. There must be some physical significance to the preference for state-change patterns lasting 30 to 60 minutes; this clearly demonstrates the need for routine long observations in order to capture and study this phenomenon properly.

3.2 Relationship to timing

The pulse period derivative for this pulsar has strong periodicities at roughly 500 and 250 days, as shown by Stairs et al. (2000) and Lyne et al. (2010). As shown by Ashton et al. (2016), the periodicity decreases roughly linearly during the timeframe covered by the dataset. In order to plot our average shape values and transition rates against a consistently-defined phase, we first convert the Lyne et al. (2010) frequency-derivative data to period derivatives $\dot{P}$, then carry out a least-squares fit with a linear slope (implying a period second derivative) as well as two harmonically-related sinusoids with linearly decreasing periods:

[TABLE]

where $P_{s}$ is the period of the cycle and $\dot{P}_{s}$ its time derivative, $t_{1}$ is a reference time, A, B, C, and D are constants, and $\phi_{2}$ is a phase offset. The fit is shown overplotted on the $\dot{P}$ data in Fig. 1. The reduced-$\chi^{2}$ of the weighted fit was 36.1, but we present the nominal uncertainties on the fit parameters in Table 5 as it is clear that the phases of the sinusoids are well-captured by this fit. We used the MJD range MJD 50870–54990 in the fit. In Fig 1 the fit is shown extended to the earlier points in the data set; if these earlier points are included in the fit, the reduced-$\chi^{2}$ is somewhat worse at 40.8, therefore we restrict the fit to the later part of the data set. For each observing epoch, we compute the phase relative to the $\sim 500$-day $\dot{P}$ cycle, taking into account the cycle period decrease. Note that phase 0 corresponds approximately to the times at which the profile is purely in the wide state. In Fig. 8 we plot the average shape and the transition rate as a function of this phase. It is clear that the pattern of the average shape is a very good match to that of the period derivative. When the profile is wide, at phase 0, the spindown rate is lowest. At later phases (near 0.3 and 0.75) in the $\dot{P}$ cycle, when the profile is narrow, the pulsar spins down faster. This reinforces the findings of Stairs et al. (2000) and Lyne et al. (2010). The arguments presented in Sec. 3.1 show conclusively that the varying widths reported in these papers were in fact due to changing fractions of time spent in each of the two emission modes. This represents a profound connection between the short- and long-term behaviours in this pulsar. Meanwhile, the rate at which transitions occur within an observation also has a $\sim$500-day periodicity, with small amounts of activity at most phases but a peak in activity around phase 0.35. While the second half of the $\dot{P}$ cycle is not quite as well-sampled, the transition rates that we do measure are generally lower. We discuss the implications of this difference in Section 4.3.

4 Discussion

4.1 Compatibility with precession

The mechanism underlying the long-term spin-down and average-profile changes in PSR B1828$-$11 is still in dispute in the literature. Free precession remains a popular hypothesis. This process involves a misalignment of the spin and angular momentum axes of the neutron star, causing the observer to have different impact angles to the spin and magnetic axes as a function of time, with the angle between these two axes remaining constant. This can naturally explain a smooth variation in pulse shape as different parts of the emission region come into view. There have been numerous attempts to model the various angles and emission-beam geometries, starting with Link & Epstein (2001), who advocated for an hourglass-shaped beam as also inferred for the geodetically precessing Hulse-Taylor binary pulsar (Weisberg & Taylor, 2002). In fact, variations on an hour-glass beam shape appear to be common to most (biaxial) precession models for this pulsar (Ashton et al., 2016, 2017), though we note that Akgün et al. (2006) derived a triaxial model with a core/cone beam shape. Intuitively, the hourglass shaped is required because the pulsar would pass from a “completely wide state” through a “completely narrow" one to a “slightly wide” one and back through the “completely narrow” state as it goes through a $\dot{P}$ cycle.

These previous beam-shape models made use of either “shape parameters” (Stairs et al., 2000) averaged over a certain range of the $\dot{P}$ cycle, or else the 10% width of the average pulse profile as determined by either a months-long average or else individual observations (Lyne et al., 2010). Both these quantities take on a continuum of values and are subject to uncertainties in the shape/width modeling for any given profile. It is therefore not surprising that a model which incorporates smooth changes in the pulse profile (that is, precession) appears in these analyses to be favoured over models allowing rapid switching between binary states (Ashton et al., 2016).

In contrast, we have demonstrated in Section 3.1 that: (1) there are abrupt, clean transitions between states in all the observations where both states were observed (modulo slight cross-contamination in the 10-s integrations on days where the switching is extremely rapid); (2) all the observations which include both states can be modeled satisfactorily as a linear combination of the two extreme states as determined from epochs with very few and completely unambiguous switches. These two properties imply that there is no difference in the overall pulse envelope from epoch to epoch. This result is completely compatible with a non-precessing pulsar for which we always see the same part of the emission region and for which the mode-switching timescales vary across the $\dot{P}$ cycle. We cannot rule out precession, but it would likely require a very small change in impact parameter in order to maintain the similar beam shape, and this fine-tuning combined with Occam’s razor leads us to strongly disfavour the precession hypothesis. It would be interesting to see the analyses of Ashton et al. (2016, 2017) redone incorporating the quantized nature of the profile switches.

4.2 Similarly switching pulsars

As shown in Lyne et al. (2010), there are now numerous pulsars which appear to undergo similar forms of sometimes abrupt switching between profile and related spin-down states. Notably, PSR B1822$-$09 also shows evidence of rapid short-term switching Lyne et al. (2010) and it is in fact known to undergo mode changes (e.g., Hermsen et al., 2017). The latter study considered 8 lengthy observing sessions showing mode changes, and did not note any differences in mode durations from epoch to epoch. The authors note, however, that all of their epochs were at times when the pulsar’s spin-down rate was similar. Other published observations of mode changing in PSR B1822$-$09 (e.g., Backus et al., 2010; Latham et al., 2012) indicate that the modes do sometimes last longer than the historical averages (Fowler et al., 1981); unfortunately the observing epochs appear to be outside the timing variation study of Lyne et al. (2010). Further long observations, triggered by changes in $|\dot{\nu}|$, could prove very informative.

Two further pulsars which have now been extensively studied are PSRs B0919+06 (Perera et al., 2015) and B1859+07 (Perera et al., 2016). Both pulsars show a “flare” profile state (Rankin et al., 2006) which may be partly related to observed changes in spin-down in PSR B0919+06 but apparently not in PSR B1859+07. Wahl et al. (2016) go so far as to attribute the flare/“swoosh” states to effects from ultra-dense binary companions. Nevertheless, the average pulse profile state does have a relationship to the spin-down in each case. Perera et al. (2015) proposed a 4-phase model for the spin-down variations in PSR B0919+06, simulating pulse arrival times and demonstrating that the frequency derivative, which necessarily must be measured across a finite span of days, can be reproduced from such abrupt changes. Ashton et al. (2016) use a similar multi-phase model when attempting to model the spin-down changes in PSR B1828$-$11 and obtained a reasonable description of the window-averaged spin-down measurements.

Whether such long-duration spin-down states with abrupt changes accurately represent the physics of this pulsar is an important question. Clearly, long-duration states with instantaneous changes in spin-down correctly describe what is seen in the long-term nulling pulsar PSR B1931+24, although there are changes in the durations of the states (Kramer et al., 2006; Young et al., 2013). However, the mode-changing data that we present in this work clearly show: (1) changes in the average shape value over time; (2) changes in the rate of mode-changing over time; (3) a relationship between both sets of changes and the phase of the $\dot{P}$ cycle. It therefore seems more plausible to infer that there are short-term changes in $\dot{P}$ corresponding to the mode changes. As $\dot{P}$ can only be measured over an observing window of weeks, its variations are smoothed out to produce the shapes seen in Fig. 1 (see Shaw et al. (2018) for an investigation of this effect). In contrast, the shape and transition rate are sampled at particular phases (Lyne et al., 2010, this work). Whether the changes in shape and transition rate are correlated at any given time – that is, tracking the smoothed $\dot{P}$ changes – or have a significant random element will have to be investigated with high-cadence observations with long tracks.

4.3 Possible chaotic behaviour

In Fig. 9 we plot the transition rate and shape data from Fig. 8 against each other, with the solid line connecting points with increasing $\dot{P}$-cycle phase. While the second half of the cycle is less well sampled, it appears that there are separate trajectories in this plane for the two halves of the $\dot{P}$ cycle, likely with some cycle-to-cycle variability. Despite this lack of coverage, it is inescapable that the point (1,0) (narrow profile, zero transitions) is reached twice in the cycle. This type of two-piece oscillation could potentially be explained by a set of nonlinear limit-cycle equations similar to the classic predator-prey equations (Lotka, 1920; Volterra, 1926). Given the broader context of the identification of chaotic behaviour in the frequency derivative of PSR B1828$-$11 (Seymour & Lorimer, 2013), we note that the shape of the plot in Fig. 9 is also reminiscent of a projected chaotic attractor, and consider whether the characteristics of the pulsar can all fit with this description.

Ashton et al. (2017) suggest that the decrease in the 500-day periodicity can be attributed to changing deformation in a precessing pulsar. However, we propose that the decrease can be accommodated within the chaotic behaviour. First, as discussed in Sec. 3.2, a linear decrease in cycle period is a poorer assumption over the full Lyne et al. (2010) data set than over the segment starting at MJD 50890. Furthermore, longer-term data will be needed to see if this trend continues. We note that PSR B0919+06 shows a decrease in its $\dot{P}$ cycle period as well (Perera et al., 2015), with a reasonably abrupt change and pattern phase shift that are not easily explained.

Seymour & Lorimer (2013) argue that there must be three variables governing the time evolution of PSR B1828$-$11. They identify the spin-down rate as one of these variables, and suggest quantities associated with the magnetic field/pulse profile for a second and the internal superfluid for the third. The average profile shape is clearly a proxy for the spin-down rate, as they are highly correlated. Here we have qualitative evidence for the transition rate between modes (which is likely related to the current and/or magnetic field properties) as a second governing variable. High-cadence data, preferably over multiple $\dot{P}$ cycles, will again be needed to explore this possibility in detail.

The two extreme profiles are well within the norm of pulsar profiles, and any single observation would simply prompt an observer to label this pulsar a mode-changer. Our observations make clear that the rate of mode changing is time-variable in a nearly repeatable fashion. Important questions are: what governs the transition rate and the average shape at any $\dot{P}$ cycle phase? If the narrow and wide extremes are both metastable states, why are there two possible “decay” trajectories from the narrow state (to average shape  0.5 and to average shape  0, with very different transition rates), and how does the system know to alternate between these two trajectories? We note that the two phases at which the profile appears purely narrow correspond to slightly different (short-term-averaged) values of $\dot{P}$. It seems inescapable that there must be yet another, as yet unknown, variable driving the system, as expected in the chaos description. Given the glitch in 2009, it seems plausible that, as Seymour & Lorimer (2013) proposed, the internal structure of the neutron star may be involved. The identification of a relationship between glitch activity and profile changes in PSR J0742$-$2822 (Keith et al., 2013) is intriguing in this respect, although the connection may not be robust (Shaw et al., in prep).

4.4 Physical mechanisms

Jones (2012) argues that the abrupt switching is due to changes in pair production, when the particle energy is pushed over a critical threshold (possibly a combination of the pair-production threshold and the work function for the neutron-star surface). This would be compatible with changes seen in ordinary moding pulsars (e.g., Wang et al., 2007). Jones (2012) further suggests that precession could cause the magnetosphere to prefer one switched state over another on the precession timescale. One possible corollary of this bias is the type of oscillation seen in PSR B1828$-$11. While we cannot definitively argue against this mechanism, we have made the point above that precession appears unlikely to be a factor for this pulsar, given the lack of evidence for continuous pulse envelope changes.

It is difficult to interpret our results concretely in terms of current theories of pulsar electrodynamics. Substantial pulse profile changes do not necessarily imply substantial changes in the global magnetic geometry. Even if the radio beam is monolithic, relatively small magnetic distortions can steer large parts of it away from the line of sight. Equally, the radio beam may be composed of multiple, overlapping sub-beams, corresponding to emission regions at various latitudes, longitudes, and altitudes in the magnetosphere (Mitra & Rankin, 2017). Relatively small changes in what field lines carry the radio-emitting conduction currents can switch some of these regions on and off, producing pulse profile changes without global magnetic changes.111This is especially true, if the emission emanates from well within the light cylinder. Near the light cylinder, the magnetic geometry is significantly modified by the Goldreich & Julian (1969) current along open field lines (and the displacement current), and it is harder to decouple current fluctuations from magnetic fluctuations. Substantial changes in the spin-down torque can also be explained without substantial changes in the global magnetic geometry. The lever arm of the torque is the light cylinder. It is therefore possible to modify the torque by redirecting conduction currents in the vicinity of the light cylinder without changing their paths very much at lower altitudes. This scenario comes with the attractive feature that it is intrinsically nonlinear, consistent with the dynamics inferred from the data, e.g., Fig. 9; current changes near the light cylinder induce magnetic changes locally (see footnote 1), which redirect the currents in turn, allowing naturally for feedback, instabilities, and multi-stable states.222Of course, just because such dynamics are allowed does not mean they occur; a specific, concrete feedback loop must first be identified self-consistently. A nonlinear version of the circuit models developed by Shibata (1991) may prove to be a profitable framework within which to analyse the possibilities. The main challenge for models of this sort is that the dynamical time-scale of the magnetosphere is the spin period, even when nonlinear effects are included. Unless some (e.g., diode-like) circuit element can be argued to function on a slow timescale, perhaps set by anomalous resistivity as in planetary magnetospheres, it is hard to get quasiperiodic changes and system memory on the timescale of years.

Long-term memory is easier to understand if it is modulated by processes in the crust, which push around the magnetic footpoints (and hence alter the global magnetic geometry), or channel conduction currents along different subsets of field lines, as discussed above. Both these scenarios would naturally lead to changes in the spin-down torque. The problem with crust-based scenarios is there is no obvious reason that they should be quasiperiodic (c.f., Huang et al., 2016). It is easy to imagine slow, stochastic drifts, analogous to polar magnetic wandering on the Earth (Macy, 1974; Wood & Hollerbach, 2015). But to get quasiperiodic dynamics, e.g. via a nonlinear limit cycle, one needs a feedback loop: the magnetospheric currents must modify the crust somehow. This is not easy to achieve, because the Lorentz forces are small compared to the elastic and hydrostatic forces that control the crust’s structure. Dynamo action is a distant possibility, but the consensus is that it is quenched by stratification (Ferrario et al., 2015, but c.f. Braithwaite (2006)).

Without trying to determine what this underlying change agent is, Timokhin (2010) proposed that relations between torque and profile may provide insights into the magnetospheric switching. Here, profile shape changes/nulling can have corresponding torque changes in the force-free magnetosphere model of pulsar emission, assuming that the magnetosphere has two or more quasi-stable states, with different closed-field regions and/or different current distributions in the open-field-line zone. An intrinsic part of that model is, however, that a large zone implies high current and energy loss. For circular beams, a larger zone implies a wider profile. Thus, according to this model, wider-profile modes spin down faster than narrow profiles. If, in pulsars like the intermittent PSR B1931+24, the narrow mode misses Earth, the spin-down during this apparent null should be less than when the pulsar is visible. That is indeed the behavior that is observed in that pulsar. In PSR B1828$-$11, however, we do see both modes. Contrary to what the Timokhin (2010) model would predict, we find the narrow profiles show the highest spin-down. The high spin-down should occur in the mode with the larger open-field-line zone. Only if the beam in the narrow mode is more elongated than in the wide mode (as seen in e.g. the fan beam producing the main pulse of PSR J1906+0746; Desvignes et al., 2013) can this larger area be realized.

In summary, therefore, there is no single scenario which addresses every serious issue adequately: nonlinear limit cycles seated in the magnetospheric suffer from a short-time-scale problem, but longer-term crustal processes are unlikely to be cyclic. The best hope for untangling the problem remains better quality polarization data, which directly probe the local (as opposed to global) magnetic geometry in the radio-emitting region.

The non-radial oscillation model (Rosen et al., 2011) may also be compatible with switching. In the case of PSR B1828$-$11, changes in the oscillation quantum number $\ell$ would provide the changes in pulse shape. Whether any of these emission models can be subject to chaotic behaviour is a subject for further theoretical investigation. Identification of additional, unobserved governing variables and subsequent full modeling of the system remains a long-term goal of these studies.

4.5 Tests for periodic magnetospheric variations

In the preferred model above, the fast profile switching and the slower $\sim$500-day $\dot{P}$ periodicity are caused by magnetospheric changes. Beyond the mode changing described above, a number of further phenomena can be used to derive the physical properties of pulsar magnetic fields. Foremost in these are measurements of subpulse drift, and the periodicity at which these drift patterns recur. Changes in drift rates are especially insightful. In van Leeuwen & Timokhin (2012), for example, these variations were used to determine relative changes in the accelerating potential drop of only 10*-3*.

We investigated whether these magnetospheric characteristics could also change over our $\sim$500-day cycle. If we could detect, in one or both of the modes, any subpulse drift or other single-pulse periodicity, then these could be tracked over the 500-day cycle, and allow for mapping of the magnetosphere and open-field-line region. We first looked if any such signals are visible in the longest, brightest sets. For these we formed uncalibrated single-pulse time series from the GBT data taken with BCPM, and we identified MJD 54127 for the narrow mode, and MJD 54202 for the wide mode, as the brightest series free of radio interference.

In each, we isolated bright consecutive single pulses which we transformed to Longitude Resolved Fluctuation Spectra (LRFS; Backer, 1970a). These Fourier transforms of consecutive pulses along each longitude can show which part (if any) of the pulse profile shows periodic behaviour such as drifting or other modulation. We took series of 2048 pulses (marked red in Fig. 3) to maximise signal-to-noise ratio while limiting the wash-out of any periodicities over longer periods. We applied the LRFS transformation using the Single Pulse Analysis toolkit (Szary & van Leeuwen, 2017), and the results are shown in Fig. 10. Compared to the large similarly derived sample in Weltevrede et al. (2006), the results lack the features expected for periodic modulation.

Given this non-detection in both modes, the prospects of using single-pulse modulation for further quantification of the $\sim$500-day magnetospheric changes are slim.

5 Conclusions

We have performed an exhaustive analysis of the variations in mode-changing seen at different observing epochs for PSR B1828$-$11. The overall pulse envelope is well-modeled by a linear combination of extreme wide and narrow profiles for each phase of the well-known $\sim$500-day $\dot{P}$ cycle. This argues against precession, which would provide us with a different view of the emission region at different $\dot{P}$-cycle phases. Histograms of state durations do not match particularly well to exponential, Poisson or Markov distributions. While there do not appear to be strong patterns within individual observations, there is some evidence for a timescale of roughly 30 to 60 minutes in data from various epochs. We compute the average mode-switching transition rate for each epoch. When plotted against the average shape value and viewed in order of increasing $\dot{P}$-cycle phase, the pattern resembles a projected double-loop chaotic attractor, matching the predictions of Seymour & Lorimer (2013) who identified chaotic behaviour in the spin-down value alone. We therefore believe we have identified this transition rate, related to the state of the magnetosphere, as a second variable in the chaotic system. We join Seymour & Lorimer (2013) in speculating that the neutron-star interior may provide the third governing variable.

Confirmation of the transition rate as the second governing variable will require considerably more data than are presented here. Lyne et al. (2010) noted that PSR B1828$-$11 frequently displayed opposite-extreme profiles in nearby observations. Our data, with observations typically lasting 4 hours, show that half-hour observations are generally inadequate to capture the true average profile (and there may well be even more transitions with timescales longer than 4 hours). Data stored in the usual pulsar timing setup, with the pulses integrated for minutes at a time, will not allow the identification of all the transitions in an observation. What is needed is multi-hour time-series observations every day or two with a telescope with the sensitivity of Parkes or the Lovell Telescope, over at least one full $\dot{P}$ cycle. Such dedicated observing time is of course hard to come by. The upcoming CHIME telescope (Bandura et al., 2014) will likely observe this pulsar on a daily basis, but it is a transit instrument operating at a frequency lower than has historically been used for this pulsar. The best hope may be to use a subarray of the Phase I SKA telescope (e.g., Karastergiou et al., 2015).

The possibility of chaotic behaviour in the mode transition rate and timing for PSR B1828$-$11 opens up the question of this pulsar’s place in the broader population. Seymour & Lorimer (2013) found patterns suggestive of but not confirming chaotic behaviour in 3 other pulsars: PSRs B1540$-$06, B1642$-$03 and B1826$-$17, namely the best-sampled pulsars presented in Lyne et al. (2010). Previous attempts to discern chaotic behaviour in pulsars, whether in timing behaviour (Harding et al., 1990) or in emission properties (DeLaney & Weatherall, 1999) did not find evidence in favour of chaos, but longer and better-sampled data sets may make a difference. As discussed above, PSR B1822$-$09 deserves further study to investigate differences in mode-duration lengths. Questions to be considered include whether the existence and timescales of these phenomena are related to the pulsar parameters such as period, age and magnetic field. In-depth studies of these objects could prove very rewarding.

Acknowledgements

IHS is supported by an NSERC Discovery Grant and by the Canadian Institute for Advanced Research. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a consolidated grant from the STFC in the UK. JvL received funding for this research from the Netherlands Organisation for Scientific Research (NWO) under project "CleanMachine" (614.001.301). DRL was supported by NSF RII Track I award number 1458952. The Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Parkes radio telescope is part of the Australia Telescope National Facility which is funded by the Australian Government for operation as a National Facility managed by CSIRO. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration, 2013). We thank Ryan Hyslop, Jennifer Riley, Raymond Lum, Cindy Tam and Sarah (Traine) Keith for their work on earlier versions of the analysis.