Independent discovery of a nulling pulsar with unusual sub-pulse drifting properties with the Murchison Widefield Array

We report the independent discovery of PSR J0027-1956 with the Murchison Widefield Array (MWA) in the ongoing Southern-sky MWA Rapid Two-meter (SMART) pulsar survey. J0027-1956 has a period of ~1.306 s, a dispersion measure (DM) of ~20.869 pc cm^-3 , and a nulling fraction of ~77%. This pulsar highlights the advantages of the survey's long dwell times (~80 min), which, when fully searched, will be sensitive to the expected population of similarly bright, intermittent pulsars with long nulls. A single-pulse analysis in the MWA's 140-170 MHz band also reveals a complex sub-pulse drifting behavior, including both rapid changes of the drift rate characteristic of mode switching pulsars, as well as a slow, consistent evolution of the drift rate within modes. In some longer drift sequences, interruptions in the otherwise smooth drift rate evolution occur preferentially at a particular phase, typically lasting a few pulses. These properties make this pulsar an ideal test bed for prevailing models of drifting behavior such as the carousel model.


INTRODUCTION
Pulsars are rapidly rotating neutron stars, the sites of some of the highest energy physical processes in the universe, owing to the very strong gravitational and magnetic fields that surround them. Discovered first via their lighthouse-like beam of radio emission (Hewish et al. 1968), which we detect as a series of regularly spaced pulses, it was quickly realised that the physics governing the relativistic plasma that generates this beam was not well understood-a state of affairs that persists to the present day (e.g. Melrose et al. 2021). The uniqueness of each pulsar's emission signature provides a wealth of information that can be used to test proposed emission mechanisms, and it has been proven again and again over the past several decades that indepth studies of individual pulsars can provide valuable clues for understanding the population as a whole. This, along with tests of general relativity (e.g. Kramer et al. 2006;Miao et al. 2021), pulsar timing arrays for gravitational wave detection (e.g. Manchester et al. 2013), pulsar braking and magnetospheric dynamics (e.g. Gao et al. 2017;Wang et al. 2020), the neutron star equation of state (e.g. Demorest et al. 2010;Antoniadis et al. 2013;Pang et al. 2021), and other pulsar science applications, motivates efforts to find new pulsars via large scale pulsar surveys, which continue to be conducted up to the present day.
One of the most recent such surveys is currently under way at the Murchison Widefield Array (MWA), a latestgeneration aperture array telescope located in Western Australia (Tingay et al. 2013). With its Voltage Capture System (VCS; Tremblay et al. 2015) and offline tied-array beamforming software (Ord et al. 2019;Mc-Sweeney et al. 2020), it has proven to be a powerful instrument for investigating a variety of pulsar phenomena such as single-pulse studies (McSweeney et al. 2017), spectral analyses (Meyers et al. 2017), studies of interstellar medium (ISM) propagation effects, and related timing applications (Bhat et al. 2018;Kaur et al. 2019). With the advent of the Phase 2 upgrade of the MWA (Wayth et al. 2018), which allows a compact configu-ration of 128 tiles (i.e. antenna elements) with short baselines ( 300 m), a pulsar survey became computationally feasible, owing to the relatively large angular size of the tied-array beam (∼ 23 ). The Southern-Sky MWA Rapid Two-metre (SMART) pulsar survey (Bhat et al. in prep) was conceived and data collection began in late 2018. At present, ∼ 70% of the data have been collected, with the initial, first-pass ("shallow") processing being limited to the first ten minutes (out of the full 80 minutes) of each observation. This strategy is beginning to pay off; the discovery of PSR J0036−1033, a lowluminosity, high Galactic latitude pulsar (b ≈ −73 • ), was reported in Swainston et al. (2021).
In this paper, we report the MWA's independent discovery of PSR J0027−1956 in the shallow pass of the SMART survey. The pulsar was originally detected in 2018 in the Green Bank Northern Celestial Cap (GBNCC) pulsar survey (Stovall et al. 2014). Unlike J0036−1033, J0027−1956 is sufficiently bright to be detected in single pulses, and was found to have a large nulling fraction. Closer inspection of J0027−1956's single pulses revealed that it belongs to the class of subpulse drifters, a phenomenon first noted in the earliest days of pulsar research (Drake & Craft 1968), and which has historically been strongly linked to magnetospheric phenomena, and thence, to the as-yet poorly understood radio emission mechanism (Ruderman & Sutherland 1975;Rankin 1986;Deshpande & Rankin 1999;McSweeney et al. 2019). "Sub-pulse" signifies discrete bursts of emission, usually much narrower than the average pulse profile, and "drifting" refers to the systematic way that sub-pulses arrive earlier or later in time in successive rotations. Sub-pulse drifting is easily detected when viewing the pulses in a stack (pulse number vs rotation phase), in which associated sets of sub-pulses form diagonal drift bands that stretch across the onpulse window. The drift rate is defined as the reciprocal of the slope of the drift bands, often expressed in units of degrees (of rotation) per pulse period ( • /P ).
Many sub-pulse drifters, e.g., B0809+74 and B0943+10, are known to have very stable drift rates (e.g. Taylor et al. 1971;Deshpande & Rankin 1999). This stability partly motivated the original carousel model, put forward by Ruderman & Sutherland (1975), which associates the sub-pulses with spatially discrete bursts of electrical discharges ("sparks") in regions of charge depletion just above the stellar surface near the magnetic poles. According to this model, the sparks are driven by E × B drift, with the observed stability of the drift rate being directly inherited from the theoretical stability of the electric and magnetic fields at the spark locations. Therefore, pulsars that exhibit anything more compli-cated than a single, stable drift rate deserve special attention, since they require modifications of, or at least extensions to, the basic carousel model. Such pulsars exist, and several extensions have been proposed over the years to account for them. For example, Gil & Sendyk (2000) suggest that a quasi-central spark can account for (non-drifting) core components in profiles, and the wellknown phenomenon of bi-drifting may be explained by presence of an inner annular gap (Qiao et al. 2004), an inner acceleration gap (Basu et al. 2020), or non-circular spark motions (Wright & Weltevrede 2017). Such extensions are typically developed to explain specific drifting behaviors that are observed in a relatively small subset of pulsars, and there still lacks a single, comprehensive theory that can describe all drifting behaviors.
Pulsars that both show complicated drifting behavior and are bright enough for single pulse analysis, are relatively rare. Even a cursory glance at the drift bands of J0027−1956 show that its drifting behavior is indeed complex and bright, making it a welcome addition to this category of pulsar. J0027−1956 is also found to contain null sequences, a phenomenon which is known to be connected to sub-pulse drifting, originally noticed by Lyne & Ashworth (1983), and further attested by recent studies of PSRs J1727−2739 (Wen et al. 2016), B1819−22 (Janagal et al. 2022), and the Vela pulsar (Wen et al. 2020).
This paper presents an analysis of the nulling and drifting behavior of J0027−1956 observed with the MWA. Section §2 describes the MWA observation in which the pulsar was originally found and the small number of archival MWA observations in which it was subsequently detected during follow up of the original candidate. Analysis of the pulsar's nulling and drifting properties follows in Section §3, which are discussed and compared to other pulsars in Section §4. We conclude with a short summary in Section §5.

OBSERVATIONS
PSR J0027−1956 was discovered in one of the observations made as part of the SMART survey (Obs ID 1226062160). A full description of the survey parameters and observational set up will be presented in an upcoming publication (Bhat et al. in prep). Many details are also described in the paper reporting the first SMART pulsar discovery (Swainston et al. 2021), for which the setup was identical, but the essential details are summarised here. The VCS delivers Nyquist-sampled dualpolarisation voltages for 128 MWA tiles (4 × 4 crossdipoles) at a rate of 100 µs and an instantaneous bandwidth of 30.72 MHz, consisting of 3072 × 10 kHz individual channels. Apart from the discovery observation, the other five observations presented in this work were taken as part of the follow-up timing campaign for J0036−1033 (Swainston et al. 2021). All six observations (summarised in Table 1) were made in the frequency range 138.88-169.6 MHz.
Each observation was individually calibrated using observations of 3C444 taken within ∼ 2 hours of the respective target observations. The calibration solutions were obtained using the Real Time System software (Mitchell et al. 2008). After calibration, the voltages were used to form a tied-array beam in the direction of the pulsar (the localization procedure is discussed in §2.1), using the bespoke beamforming software described in Ord et al. (2019). The MWA is located in a very radio quiet location, and no significant radio frequency interference was found in the beamformed data. After determining the best dispersion measure (DM) and period (see §2.2), dedispersion and folding were performed with DSPSR (van Straten & Bailes 2011) and PSRCHIVE (Hotan et al. 2004) to form single pulse archives, from which Stokes I pulse stacks were created. Profiles of each detection, as well as an integrated profile incorporating the data from all four observations, is shown in Fig 1. Subsequent analysis of the pulse stacks is described in §3.
Although this pulsar was detected in the shallow pass of SMART, it was subsequently noted to have been reported in 2018 as a candidate in the GBNCC pulsar survey (Stovall et al. 2014). Later re-processing of data from the 70-cm Parkes Southern Pulsar Survey (Manchester et al. 1996) also revealed a weak detection that was not previously reported.
Follow-up of the initial MWA detection with the full 80 minutes, as well as in other archival MWA observations, revealed that the pulsar has a very large nulling fraction (∼ 77%) 1 , with some nulling sequences (excluding occasional intermittent pulses) potentially exceeding 20 min (i.e. the length of the first observation on MJD 59002, which contained only a few intermittent pulses, as shown in Fig. 2). A full analysis of the nulling properties is presented in Section §3.1 below. Despite the fact that J0027−1956 was serendipitously discovered in the first ten minutes, it may well be the harbinger for a population of relatively bright, intermittent pulsars which the SMART survey, with its 80-minute dwell times, is well-suited to detect. A similar strategy is employed in the Low Frequency Array (LOFAR) Tied-Array All-Sky Survey (LOTAAS; Sanidas et al. 2019) and the low latitude (|b| ≤ 5 • ) segment of the High Time Resolution Universe survey (HTRU; Keith et al. 2010), which use dwell times of approximately one hour.

Localization
The discovery observation (Obs ID 1226062160) was taken when the MWA was in its compact configuration (see Wayth et al. 2018), giving an effective tied-array beam size of ∼ 23 at the central frequency of 155 MHz.
As per the SMART survey design, these beams are tiled across the field of view, each of which is searched for periodic candidates. The compact configuration includes many redundant baselines, and the tied-array beam for the compact configuration is complex. Consequently, the pulsar was detected in multiple adjacent beams, including a boresight beam that yielded the highest signalto-noise, and several nearby grating lobes. The ability to exploit grating lobe detections for fast-tracking candidate confirmation and localization will be discussed in the survey description paper (Bhat et al. in prep).
The other five observations were taken when the MWA was in its extended configuration, giving a tied-array beam size of ∼ 3 at 155 MHz. Localization within the four observations in which the pulsar was bright enough followed the identical procedure as used in Swainston et al. (2021), i.e. gridding the candidate positions with overlapping tied-array beams and using signal-to-noise and the estimated beam shape to produce a likelihood map of the source position. It is known that the ionosphere is capable of shifting the apparent position of sources by an appreciable fraction (∼ 10-30%) of the tied-array beam, and that these offsets are not necessarily corrected by the calibration process when the calibration observations are separated from the target observation, both spatially and temporally (Swainston et al. 2022). Nevertheless, the localization afforded by the extended array observations was sufficiently precise (a few arcminutes) to warrant follow-up with a highresolution imaging telescope. Accordingly, images made with Band 3 (300-500 MHz) data from the upgraded Giant Metrewave Radio Telescope (uGMRT; Gupta et al. 2017), obtained under Director's Discretionary Time (DDTC185), enabled us to confirm the pulsar's position to an accuracy of ∼ 1 . The right ascension and declination of the final best position, consistent with both the uGMRT imaging and the MWA extended configuration observations, are 00h26m36s and −19 • 55'59".

Determination of Period and Dispersion Measure
The limited number of available detections make obtaining a timing solution impossible at this early stage. Further follow up observations with the uGMRT and the Parkes 64-m radio telescope (also known as Murriyang), as well as continued processing of archival MWA observations, are currently underway. These observations, along with the full timing solution that they are expected to yield, will be reported in a future publication.
Lacking a timing solution, we used PSRCHIVE's pdmp routine to determine the best period (P ) and DM for each observation independently. This routine performs a brute-force grid search in both period and DM pa-rameter space, and returns the period and DM that yields the highest profile signal-to-noise. The weighted averages of both period and DM (respectively) over the four observations were calculated, with the uncertainties reported by pdmp being used to weight the individual measurements. The values thus derived were P = 1.306150±0.000005 s and DM = 20.869±0.005 pc cm −3 . This simple estimate does not take the period derivative into account (in fact, it implicitly assumesṖ = 0), and we note that a sufficiently largeṖ would cause the period to drift significantly over the time period spanned by our four MWA observations (∼ 1.5 yr). However, the pulse stacks (e.g. Fig. 4) folded on the above period showed no detectable slope, which we estimate would be discernible if the pulse window changed by more than 20 ms over the course of the observation. Given this, we place a conservative upper limit of |Ṗ | ≤ 10 −13 s s −1 .
Any error in the folding period also naturally translates into a systematic error in the sub-pulse drifting analysis presented below, e.g. the drift rate. However, the fractional uncertainty of the period is much smaller than those of the sub-pulse drifting analysis, which are dominated by the stochastic and bursty nature of the sub-pulses themselves. Hence, in the following analysis, we neglect this systematic error, which is not expected to affect the main results.

Nulling properties
The emission from J0027−1956 appears either as long ( 10 min) burst sequences exhibiting phase-modulated sub-pulse drifting, or as short ( 20 s) bursts interspersed throughout otherwise long (up to 25 min) nulling periods. Both types of emission are visible in Fig. 2, in which the peak flux density within the pulse window of each pulse is plotted as a function of time for all four MWA observations.
To estimate the nulling fraction, we first identified which pulses contained sub-pulses within the pulse window (this is a semi-automated procedure described in §3.2). Any pulse without a sub-pulse is then counted as a null, even if it is nested within a burst sequence. Under this definition, we report a nulling fraction of ∼ 77%, calculated using all six MWA observations (refer to Table 1). Note that if the two observations taken on MJD 59002 are excluded, the nulling fraction would be ∼ 72%. This working definition of a null may misclassify pulses with broad emission spread out over the pulse window, even if its integrated flux density is statistically significant. For comparison, we present the pulse energy histograms of both on-and off-pulse regions of the pulse stacks for each of the six observations in Fig. 3. The similarity between the on-and off pulse distributions for observations 1275094456 and 1275172216 indicate that any broad emission that is present falls well below the MWA detection threshold. Despite this, the uncertainty in our estimate of the nulling fraction will still be dominated by the lack of a statistically significant number of burst sequences.
The error on this estimate, however, will be dominated by the small number of burst/null sequences captured in these six observations. Therefore, the true nulling fraction may ultimately prove to be much lower or higher than the ∼ 77% quoted here. Ongoing follow-up observations (and a more complete search through archival MWA data) will allow us to place much stronger constraints on the nulling fraction, and will be reported in a future publication.

Drifting behavior
During burst sequences, the single pulses form distinct phase-modulated drift bands that are always oriented such that sub-pulses associated with the same drift band appear at earlier rotation phases over time ('positive drifting' in the classification scheme of Basu et al. 2019). The drift rate, however, is highly variable, as can be seen in the example pulse stacks in Fig 4. Sometimes the drift rate appears to change abruptly, reminiscent of the drift mode changes observed in several other pulsars (e.g. Kloumann & Rankin 2010;Rankin et al. 2013;Joshi et al. 2017;McSweeney et al. 2019). However, apart from these abrupt changes, the drift rate is seen to evolve gradually, as does the "vertical spacing" between them, denoted by P 3 (see, e.g., the sequences marked "A1" in Fig 4). The most dramatic examples of this slow evolution show a difference of more than a factor of two between the value of P 3 at the beginning and end of a burst sequence.
The manner in which the drift rate and P 3 evolve make the clean division of burst sequences into distinct drift modes difficult. Nevertheless, the most economical description of the observed drifting behavior is that it consists primarily of two modes, which we label A and B. We further divide these into two subcategories (A1/A2 and B1/B2) which depend on qualitative properties of their appearance and context. This categorisation scheme is described below.

Mode classification
Mode A is characterised by slowly evolving P 3 values that can range from ∼ 20 to ∼ 55 P , and drift rates in the range −0.3 to −0.7 • /P . Most commonly, P 3 is seen to increase over time, as the examples in Fig. 4 show. However, there is at least one clear instance in the (MJD 59094 observation where P 3 decreases from approximately 35 to 24 P over the course of 170 rotations.
Mode A sequences can last anywhere from a few tens to a few hundreds of pulses. For longer sequences, Mode A is easily identifiable from the relatively large P 3 , but for shorter sequences, only one or two drift bands are visible, and therefore P 3 cannot be measured directly. In these cases, a positive identification of Mode A is based on the drift rate, which can be measured even for single drift bands.
The distinction between the subcategories A1 and A2 is based primarily on the level of organisation of the drifting behavior. Sequences whose drift bands follow  smooth (albeit curved) tracks with little deviation are classified as A1. On the other hand, A2 sequences show frequent interruptions, both in terms of short-duration nulls as well as rapid, but temporary, deviations in the drift rate. These interruptions are observed to last between ∼ 5-10 pulses, but afterwards the drift bands reappear at the rotation phase that one would expect if the interruptions had not occurred. These are further discussed in §3.2.3.
Mode B sequences are relatively short (up to 30 pulsar rotations in duration), and have small P 3 values (∼ 10 P ) and large drift rates (∼ −1.2 • /P ). Longer sequences ( 10 rotations) are found to occur most commonly nearby to Mode A sequences; these are classified as B1. Shorter sequences ( 10 rotations) may be found anywhere, even in the midst of long, otherwise uninterrupted null sequences. As with Mode A sequences, those sequences which are too short in duration to admit more than one or two drift bands may still be positively identified by their drift rate.

Drift band modelling
The pulse stacks in Fig 4 clearly show that the slow evolution of the drift rate is nevertheless "fast" enough that even individual drift bands (which can span in excess of 50 pulses) can show significant curvature. Following the lead of Lyne & Ashworth (1983), we therefore modelled the drift bands with an empirical function that assumes an exponential decay rate for the drift rate, D, following a null: where D 0 is the difference between the asymptotic drift rate, D f , and the drift rate at the onset of the drift sequence; p is the number of pulses since that onset; and τ r is the drift rate relaxation time (in units of the rotation period). The fitting procedure is carried out as follows. The beginning and ending pulses of a drift sequence were manually identified by visual inspection of the pulse stack. Sub-pulses were then identified by first smoothing the pulses with a Gaussian kernel of width ∼ 3.6 ms (i.e. 1 • of pulsar rotation, the approximate width of a subpulse), and identifying peaks above a certain flux density threshold. The threshold is chosen to be the minimum pixel value for which no sub-pulses are identified in the off-pulse region throughout the entire pulse stack. Each sub-pulse in the drift sequence is then assigned a drift band number, and fitted to the following functional form of the sub-pulse phases (using SciPy's curve fit method). This function is obtained by considering the drift rate as the rate of change of sub-pulse phase with pulse number, and integrating Eq. (1): where ϕ is the phase of a sub-pulse; ϕ 0 is an initial reference phase; d is the (integer) drift band number; and P 2 is the longitudinal spacing between successive drift bands. In total, this model has five free parameters, D 0 , D f , τ r , ϕ 0 , and P 2 , of which the expression ϕ 0 +P 2 d defines the pulse phase at p = 0. During the fitting, no restriction is placed on the sign of τ r . A positive value would mean that the drift rate eventually asymptotes to D f , whereas a negative value would mean that the drift rate is increasing exponentially from an initial drift rate (i.e. D → D f as p → −∞). Unless the relaxation time is much less than the duration of the drift sequence, τ r is highly covariant with the other parameters, and the model is not capable of distinguishing between the two scenarios with any confidence. Given this, the primary utility of the empirical model is its ability to predict other measurable quantities such as P 2 , P 3 , and the drift rate, which are less sensitive to the precise functional form used.

P2, P3, and drift rate
The drift band fitting model defined in Eqs. (1) and (2) can be used to extract information about P 2 , P 3 , and the drift rate, D. P 2 , which is assumed to be constant everywhere throughout a drift sequence, is fit independently for each drift sequence. P 3 and the drift rate are defined continuously over all pulse numbers (even fractional ones) by virtue of the fact that the drift bands are modelled as continuous functions on the pulse stack, although we caution that these quantities are only actually measured within the pulse window, at intervals of a near-whole number multiple of the rotation period. Eq. (1) gives the drift rate, and one can then derive P 3 = P 2 /|D| directly. The assumption of constant P 2 (for a given drift sequence) is justified by virtue of the goodness-of-fit of the models to the drift bands. Fig. 5 shows the result of the drift band fitting to the MJD 58991 observation, alongside the values of P 2 , P 3 , and the drift rate predicted by the model. One can immediately note several interesting features. First, the modelling bears out the observation that P 3 can indeed change gradually over the course of the mode A drift sequences by a large amount, even by more than a factor of two (from ∼ 20 to 40 P in one case). Second, P 2 appears to be generally smaller for mode B sequences (∼ 12 • ) than mode A (∼ 14 • ). Whether or not this is true generally, or whether a single value of P 2 could be used to model all drift sequences simultaneously has not been attempted here, but in any case would require a larger sample of drift sequences to test robustly.
Third, the nature of the disorganization of A2 sequences can now be more closely examined. Fig. 6 shows such an example A2 sequence to which the above model has been fitted. Although the model fits well overall, non-random deviations from the fit can be seen in several drift bands. These deviations seem to consist of an initial rapid increase of the drift rate (i.e. the drift bands appear to veer to the left), followed by a short null (only a few pulses long), after which the pulses reappear on the trailing side (i.e. the right hand side) of the fitted model before rapidly approaching the fitted model once again. The temporarily higher drift rate appears visually similar to the drift rate observed during mode B, suggesting that these interruptions are actually brief excursions into mode B. Curiously, all such deviations (across all observations) appear at about zero phase, which approximately corresponds to the peak of the profile.

Fluctuation spectra
Fluctuation spectra are most useful when the drift bands are regularly spaced (i.e. when P 3 is constant). The fact that this is clearly not the case for J0027−1956 means that features in the fluctuation spectra will be spread out over a range of spectral bins. Fig. 7 shows the longitude-resolved fluctuation spectrum (LRFS) and the two-dimensional fluctuation spectrum (2DFS) for the drift sequence shown in Fig. 6. Despite the changing P 3 and the relative disorganization of the drift rate, the LRFS and 2DFS show surprisingly high-Q features at P 3 ≈ 33±4 P and P 2 ≈ 16 • ±2 • , which are consistent with the fitted model values of P 3 changing from ∼ 50 P to 27 P throughout the sequence, and P 2 = 13.8 • ± 0.3 • . The relative lack of power in the bottom panel of the LRFS near phase ∼ 2 • is likely related to the interrup- Figure 6. A sequence of mode A2 from the MJD 58434 observation, with fitted model drift band overlaid in white. The "interruptions" from the otherwise good fit can be seen most clearly around pulse numbers ∼ 450, ∼ 490, and ∼ 620.
tions in the organization of the drift bands that occur at that phase, discussed in the previous section. Interestingly, the presence of extra power in the LRFS in the leading component between 0.05-0.1 cycles/P is suggestive of a connection to Mode B, which (with its P 3 ≈ 10 P ) would be expected to contribute power in the LRFS near this range.

DISCUSSION
PSR J0027−1956 is a clear example of a pulsar that exhibits sub-pulse drifting; its drifting behavior is qualitatively similar to other pulsars of the same class. The slope of the drift bands always has the same sign (later sub-pulses arrive at earlier rotation phases than their predecessors), and the drift bands are almost universally connected across the entire pulse window.
These basic facts, when interpreted in the context of the hollow cone/carousel model (Ruderman & Sutherland 1975;Rankin 1986;Deshpande & Rankin 1999), suggest that J0027−1956 is a "conal single" pulsar, where the line of sight cuts tangentially through the cone of emission that plays host to a carousel of discrete beams rotating around the star's magnetic axis. The profiles (see Fig. 1) show a barely-resolved double peak structure, which suggests that the line of sight is just on the poleward side of the emission cone at this frequency.
The most remarkable feature of J0027−1956's drifting behavior is how the drift rate evolves over time, exhibiting both rapid changes (indicating the presence of  Fig. 6. Bottom: A zoomed section of the 2DFS of the same pulse stack section. The feature at ∼ 0.03 cycles/P corresponds to P3 ≈ 33 ± 4 P , consistent with the average P3 of this sequence determined from the pulse stack directly. The more diffuse power in the range 0.05-0.1 cycles/P in the leading component (and the corresponding feature visible in the 2DFS) are indicative of the relative disorganization of the drift bands that occurs there. multiple, distinct drift modes) as well as gradual, intramode changes. The presence of multiple drift modes alone is known to occur in several other pulsars, e.g. B0031−07 (Huguenin et al. 1970;Joshi & Vivekanand 2000) and B1944+17 (Deich et al. 1986;Kloumann & Rankin 2010), but unlike J0027−1956, the drift rates of these pulsars' modes are usually stable enough that the modes can be uniquely characterised by their P 3 values.
Because of this inherent stability, analyses of how the drift rate changes within drift modes (or if it changes at all) are relatively rare, but several types of intramode evolution of the drift rate have been historically identified. B0809+74, for example, temporarily adopts a slightly lower drift rate immediately following a null, which then relaxes slowly back to the usual, stable rate after a few tens of pulses. Lyne & Ashworth (1983) demonstrated that this variability is related to the duration of the null sequence preceding it, and that the drift band phase is "remembered" across the null. B0826−34, in contrast, has a drift rate that varies slowly but pseudo randomly at all times, making the identification of discrete drift modes impossible (e.g. Esamdin et al. 2012). More recently, McSweeney et al. (2017) showed that B0031−07, which has long been known to exhibit three stable modes, has a measurable slow evolution of the drift rate within individual drift sequences. The properties of these and a few other, similar pulsars are presented with J0027−1956 in Table 2.
Of the pulsars listed in the table, J0027−1956 is arguably most similar to B0031−07, with which it shares (1) the presence of more than one drift mode, (2) the fact that P 3 changes slowly throughout individual drift sequences (McSweeney et al. 2017), and (3) a relaxation time (if such a model is even applicable) that is longer than the typical duration of the drift sequences themselves. This last point is what renders the particular drift band model used in this work unable to measure the relaxation time reliably; in this respect it is not significantly better than the "quadratic" model used in Mc-Sweeney et al. (2017), which is equivalent to keeping only the lowest order terms in the Taylor expansion of Eq. (2). Notably, Joshi & Vivekanand (2000) report that in B0031−07, drift phase is indeed remembered across nulls, as first reported for B0809+74 and B0818−13 by Lyne & Ashworth (1983), and a visual inspection of the J0027−1956 pulse stacks suggest that a similar phase memory is also present here, although a more careful analysis is required before this can be verified.
In general, the connection between nulls and drift rate is well established, but the nature of this relationship varies from pulsar to pulsar (see Table 2 for references). The exponentially decaying drift rate reported for B0809+74 has a relaxation time of between 10 and 20 pulses, while B0818−13's drift rate varies like a damped oscillator on a similar timescale. In contrast, fluctuation spectral analyses of B0943+10 reveal that its drift rate evolves very slowly after a null, only settling to its stable value after several thousand pulses. B0826−34's drift rate evolves on a faster time scale, but in a pseudorandom manner, not unlike B2016+28, whose drift rate can vary on a time scale of a few tens of pulses.
All of these complex drifting behaviors present a natural challenge to the carousel model, whose most basic formulation predicts only a single stable drifting mode, due to the expected stability of the underlying elec-tric and magnetic fields that generate the carousel motion. Nevertheless, many of the above effects can be successfully explained with the carousel model by invoking aliasing effects that arise due to beating between the stellar rotation rate and the carousel's rotation. In the presence of aliasing, multiple drift modes can be explained if the number of sparks in the carousel changes, even if the carousel speed itself doesn't change. This idea has been invoked to explain the multiple drift rates of B1918+19 (Rankin et al. 2013), B1819−22 Janagal et al. 2022), and B0031−07 (Mc-Sweeney et al. 2019). Aliasing can also "amplify" the variability of the drift rate if the ratio of stellar rotation and carousel rotation is very close to a ratio of small integers. That is, a small change in the magnetospheric conditions (e.g. the local electromagnetic configuration) near the stellar surface (in the gap, where the sparks reside) can appear to the observer as a relatively large change in the drift rate. van Leeuwen & Timokhin (2012) elucidated the connection between the accelerating potential within the polar gap and the carousel behavior, and found that the variability of B0826−34's drift rate can be explained in this way. In contrast, van Leeuwen et al. (2003) argued that aliasing could not be present during the post-null changing drift rate for B0809+74. Whether or not aliasing is included, models for changing drift rates are often discussed in terms of the local magnetospheric conditions changing over the relevant time scales (e.g. van Leeuwen et al. 2003;Yuen 2019).
One must question whether or not these varied behaviors are all manifestations of the same underlying phenomenon, operating on different timescales. That is, the general rule that governs all these pulsars might be that the drift rate varies in a systematic way following (and sometimes preceding) a null sequence, and that given ample time before another interruption, all would settle into a stable drift rate. This is difficult to test for pulsars like J0027−1956 and B0031−07 unless a sufficiently long, uninterrupted drift sequence is observed with a sufficiently short relaxation time. Among the observations of J0027−1956, only one mode A1 drift sequence was found whose fitted relaxation time of τ r ≈ 31 P was significantly shorter than the drift sequence itself (285 pulses). This sequence starts with P 3 ≈ 12 P , and after stabilizing, continues on for well over 100 pulses with P 3 ≈ 40 P .
Modelling the drift sequences in this way opens up new avenues for exploring and quantifying the relationships between the drift rate, the drift decay rate, the duration of the preceding drift (or null) sequence, the possible phase connection between consecutive drift modes (or across nulls), and others. For example, in many instances the drift bands appear to connect smoothly across mode switches between modes A and B, which, if true, has implications for the mechanisms by which mode transitions occur in the context of the carousel model. Related to this is the possibility that the drift rate "interruptions" seen in mode A2 sequences are actually miniature excursions into mode B. If so, why these excursions occur preferentially within a particular phase range also has implications for the carousel model. The verification and exploration of such claims require much larger data sets than are currently available, and will be explored in a follow-up paper with an enlarged data set comprising more archival MWA observations as well as recently taken uGMRT observations of this pulsar.

CONCLUSIONS AND SUMMARY
We have reported the independent discovery and first MWA observations of PSR J0027−1956 in the SMART pulsar survey. It is bright enough to see in single pulses at MWA frequencies, but contains long null sequences which in some cases can apparently exceed 20 minutes. Such pulsars are inherently difficult to detect because of their intermittent nature, but are expected to be discovered in greater numbers in surveys like SMART and LOTAAS, which employ long ≥ 1 hr dwell times.
Its slow period and intermittent nature mean that long duration observations are required to collect a sufficiently large number of complete burst sequences to undertake robust statistical analysis on their properties. Follow-up observations are currently under way at the uGMRT and at Murriyang, which will allow us to explore the drifting behavior of J0027−1956 in much greater depth, as well as obtain a timing solution, refine the nulling fraction estimate, and perform a polarimetric analysis of its single pulses.
However, even with the relatively few observations presently at our disposal, J0027−1956 clearly exhibits many interesting single pulse behaviors. Among these, the most prominent is the manner in which the drift rate changes slowly throughout drift sequences (a property observed in only a small number of other sub-pulse drifting pulsars), as well as the interplay between the two drift mode classes identified here (modes A and B) and the intervening null sequences. These properties make the interpretation of this pulsar's drifting behavior particularly interesting in the context of the carousel model. In-depth investigations of this kind will thus play an important role in the quest to uncover the intricacies of pulsar emission physics, an outstanding problem in pulsar astronomy. This scientific work makes use of the Murchison Radioastronomy Observatory, operated by CSIRO. We acknowledge the Wajarri Yamatji people as the traditional owners of the Observatory site. Support for the operation of the MWA is provided by the Australian Government (NCRIS), under a contract to Curtin University administered by Astronomy Australia Limited. We acknowledge the Pawsey Supercomputing Centre which is supported by the Western Australian and Australian Governments. This work was supported by resources awarded under Astronomy Australia Ltd's AS-TAC merit allocation scheme on the OzSTAR national facility at the Swinburne University of Technology. The OzSTAR program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government. The GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research, India. Software support resources awarded under the Astronomy Data and Computing Services (ADACS) Merit Allocation Program. ADACS is funded from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government and managed by Astronomy Australia Limited (AAL). DK was supported by NSF grant AST-1816492. We would also like to thank Alex McEwen for making the original discovery of this pulsar in the GBNCC survey. RMS acknowledges support through Australian Research Council Future Fellowship FT190100155. We thank the anonymous referee for the many detailed suggestions which improved this paper.