An Unusual Pulse Shape Change Event in PSR J1713+0747 Observed with the Green Bank Telescope and CHIME

The millisecond pulsar J1713+0747 underwent a sudden and significant pulse shape change between 2021 April 16 and 17 (MJDs 59320 and 59321). Subsequently, the pulse shape gradually recovered over the course of several months. We report the results of continued multifrequency radio observations of the pulsar made using the Canadian Hydrogen Intensity Mapping Experiment and the 100 m Green Bank Telescope in a 3 yr period encompassing the shape change event, between 2020 February and 2023 February. As of 2023 February, the pulse shape had returned to a state similar to that seen before the event, but with measurable changes remaining. The amplitude of the shape change and the accompanying time-of-arrival residuals display a strong nonmonotonic dependence on radio frequency, demonstrating that the event is neither a glitch (the effects of which should be independent of radio frequency, ν) nor a change in dispersion measure alone (which would produce a delay proportional to ν −2). However, it does bear some resemblance to the two previous “chromatic timing events” observed in J1713+0747, as well as to a similar event observed in PSR J1643−1224 in 2015.


INTRODUCTION
PSR J1713+0747 is a bright millisecond pulsar (MSP) that serves as a sensitive component of pulsar timing arrays (PTAs), projects that aim to detect gravitational waves at nanohertz frequencies through their influence on the arrival times of pulses from a large number of MSPs over many years.Its brightness (approximately 8.3 mJy at L band; Spiewak et al. 2022a) and the sharp features in its pulse profile (effective width 0.54 ms at L-band; based on data from Alam et al. 2021a) make possible particularly high-precision time-of-arrival (TOA) measurements.As a result, J1713+0747 is currently observed by all member collaborations of the International Pulsar Timing Array (IPTA; Perera et al. 2019), namely, the North American Nanohertz Observatory for Gravitational Waves (NANOGrav; Alam et al. 2021a;Alam et al. 2021b), the European Pulsar Timing Array (EPTA; Desvignes et al. 2016;Chen et al. 2021), the Parkes Pulsar Timing Array (PPTA; Kerr et al. 2020;Reardon et al. 2021), and the Indian Pulsar Timing Array (InPTA; Joshi et al. 2018;Nobleson et al. 2022), as well as by participating scientists associated with the Five-hundredmeter Aperture Spherical Telescope (FAST) in China (Nan & Li 2013;Hobbs et al. 2019) and the MeerKAT telescope in South Africa (Spiewak et al. 2022b).In 2013, Dolch et al. (2014) undertook a 24-hour global observing campaign, using nine IPTA telescopes to assess the precision timing capabilities of PSR J1713+0747.They found that TOA measurement precision was ultimately limited by pulse jitter to a level of 27 µs for single pulses, which would imply a TOA error of 6.3 ns for the entire 24-hour observation.In addition to its use by PTAs, the high timing precision achievable with J1713+0747 and its well-characterized relativistic orbit make it an important pulsar for tests of general relativity (e.g.Zhu et al. 2019).
In April 2021, a large pulse shape change, noticeable by eye in comparisons of profiles, was detected in J1713+0737 (Xu et al. 2021).This change was surprising because the mean pulse profiles of pulsars are generally observed to be constant over long periods of time (e.g., Cordes 2013).They presented profiles from before and after the shape change, observed with the Effelsberg and Nançay radio telescopes at C band (4.9-5.1 GHz) and L band (1228( -1740 MHz) MHz), respectively, as well as with the Kunming 40-meter telescope at S band (2150-2400 MHz) and with FAST at L band (1050-1450 MHz).The FAST profiles included polarization information, indicating that the linearly polarized emission from the pulsar also changed as a result of the event.Additionally, Xu et al. (2021) noted that the shape change appeared to be accompanied by a change in the pulsar's dispersion measure (DM) of approximately −4.3 × 10 −3 pc cm −3 , which they measured using the FAST L-band data.The existence and timing of the event were subsequently confirmed using other telescopes, including the Giant Metrewave Radio Telescope (GMRT; Singha et al. 2021), the Green Bank Telescope (GBT; Lam 2021), and the Canadian Hydrogen Intensity Mapping Experiment (CHIME; Meyers & CHIME/Pulsar Collaboration 2021).Since CHIME observes the pulsar with daily cadence, Meyers & CHIME/Pulsar Collaboration (2021) were able to constrain the time of the event's onset to a 24-hour period, between MJD 59320 and MJD 59321 (April 16 and 17, 2021).In the months after the shape change, the frequency-dependent pulse shape and arrival times have been observed to recover toward typical values seen prior to 2021 (Jennings et al. 2022).
Importantly, the 2021 shape change event cannot be explained as a change in DM alone.A change in DM would not produce the complex changes in pulse shape that are observed, and the frequency dependence of the TOA residuals deviates significantly from the ν −2 scaling expected from dispersion.
On the other hand, this new shape change event in J1713+0747 is not entirely without precedent.PTA observations of J1713+0747 have shown two unusual chromatic (i.e., radio-frequency dependent) timing events in the last 15 years, although both are smaller in amplitude than the most recent event by approximately an order of magnitude.The first of these events, which took place around MJD 54750 (October 11, 2008), was identified by Demorest et al. (2013).The second event took place around MJD 57510 (May 2, 2016), and was identified by Lam et al. (2018).Both events show an abrupt onset, lasting no more than a few days, after which pulses appear to arrive earlier.This is followed by a gradual, approximately exponential recovery taking place over several months.Since the TOA measurements change more at lower frequencies, each event produces a change in apparent DM.For the first event, this is about −6 × 10 −4 pc cm −3 , while for the second event, it is about −4.3 × 10 −4 pc cm −3 .Lam et al. (2018) attribute the events to lensing of the radio emission by some structure in the ionized ISM.Follow-up by Lin et al. (2021), also using the NANOGrav 12.5-year data, showed that the data are consistent with a model in which the lensing is produced by a region of lower electron density with the geometry of a folded sheet.Shape changes associated with the first two chromatic timing events are not apparent by eye, but Goncharov et al. (2021) and Lin et al. (2021) have found that they do occur at a low level.
Another form of transient pulse shape change was observed in the Crab pulsar, PSR B0531+21, in 1997, where it was attributed to refractive lensing of the pulsar's radio emission by a sharp feature in the density distribution of plasma at the edge of the Crab nebula (Backer et al. 2000;Graham Smith et al. 2011).Given this context, it is apparent that there are many possible explanations for the event considered here.Several of these possibilities will be considered in more detail in Section 4 below.
In what follows, we report on recent observations of J1713+0747 made by NANOGrav and the CHIME/Pulsar Collaboration using the GBT and CHIME, both before and after the shape change event, and discuss some possible astrophysical scenarios which may have given rise to the event.In Section 2, we describe the cadence and frequency coverage of the observations, the receivers used to make them, and the analysis techniques used to reduce them.
In Section 3, we analyze the main features present in the data.Then, in Section 4, we consider several possible astrophysical interpretations for the observed phenomena.Finally, in Section 5, we summarize our results and conclude.

Observation
The observations considered here were made using the GBT and CHIME telescopes.The GBT observations were made using two receivers: the 800 MHz band of the PF1 prime focus receiver (PF1/800), which has a bandwidth of 200 MHz centered on 820 MHz, and the L-band Gregorian focus receiver, which has 800 MHz of bandwidth centered on 1500 MHz.The observations were coherently dedispersed and folded in real time using the VEGAS backend (Roshi et al. 2011;Prestage et al. 2015) in its pulsar mode.Pulse profiles with 2048 phase bins were produced in initial channels 1.5625 MHz wide, with 512 such channels across the band at 1500 MHz and 128 at 820 MHz.GBT observations at both 820 and 1500 MHz were made roughly monthly as part of the main NANOGrav observing program, and typically have integration times of about 30 minutes.Additional 1500 MHz observations of J1713+0747 were formerly made as part of the NANOGrav high-cadence observing program; unfortunately, these higher-cadence observations ended in March 2021, approximately a month before the shape change event began, so the average cadence of 1500 MHz observations is lower after the onset of the shape change than before it.For a short (approximately 1 minute) period before each GBT observation, an artificial 25 Hz pulsed noise source was injected into the signal path and recorded.These noise diode observations were later used to calibrate the differential gain and phase between the two hands of polarization of the receiver (cf.NANOGrav Collaboration et al. 2015, section 3.1).
The CHIME observations were made using the telescope's 256-antenna receiver system, which has 400 MHz of bandwidth centered on 600 MHz (CHIME Collaboration et al. 2022).Since CHIME consists of four stationary cylindrical reflectors oriented in a north-south direction, it observes a stripe along the meridian 120 • long and 1.3 • -2.5 • wide at any given moment, and cannot be pointed.Instead, the CHIME/Pulsar backend uses digitally synthesized beams to observe as many as ten pulsars simultaneously as they transit the main beam (CHIME/Pulsar Collaboration et al. 2021).This makes it possible to observe nearly every pulsar visible to CHIME once every day, but only for a limited time, which for J1713+0747 is typically between 13 and 15 minutes.In practice, J1713+0747 was observed every day as it transited, except for the relatively few days when the telescope was shut down, either for maintenance or because of extreme high temperatures (more than 45 • C).This gives the CHIME observations much higher cadence than even the "high-cadence" GBT observations.Like the GBT observations, the CHIME observations were also coherently dedispersed and folded in real time, to produce pulse profiles with 1024 phase bins in each of 1024 initial channels.However, due to the very different nature of observing with CHIME, noise diode observations were not available, and the polarization calibration performed for the GBT data could not be carried out (cf.CHIME/Pulsar Collaboration et al. 2021, section 4.4).For this reason, we do not attempt to interpret the polarization properties of the CHIME data below, but consider only the total intensity data.
A summary of all observations used in this paper is given in Table 1.Below, we consider primarily those observations made between February 21, 2020 and February 1, 2023 (MJDs 58900-59976); observations from before this period were used only in constructing templates.GBT made a total of 29 observations with the 820 MHz receiver and 70 observations with the 1500 MHz receiver during this period.Of these, 16 of the 820 MHz observations and 21 of the 1500 MHz observations took place after the shape change event.In the same period, CHIME made 884 observations of the pulsar, with 541 taking place after the shape change event.Due to the recentness of the event, none of the observations used in this paper are part of NANOGrav's published 12.5-year data set or the 15-year data set currently being prepared1 .However, all of them will be included in a future NANOGrav data set.Additionally, the full set of data used in this paper, including the frequency-resolved profiles derived from each observation, is available in an accompanying Zenodo dataset 2 .

Post-processing and alignment
To prepare them for analysis, the GBT data were processed using NANOGrav's nanopipe software (Demorest 2018), which excised pre-determined bands containing significant radio-frequency interference (RFI) or receiver resonances, and performed a basic polarization calibration procedure.No such pre-determined bands for RFI excision were available for CHIME, nor were the CHIME observations accompanied by dedicated calibration scans, so the CHIME data were not polarization calibrated, and all RFI removed was identified manually at a later step.All the data were then dedispersed at a constant reference DM of 15.9638 pc cm −3 and aligned using a pulsar ephemeris fit to the NANOGrav 15-year data (Agazie et al. 2023).The exact parameter values used can be found in the Zenodo dataset accompanying this paper.Finally, they were post-processed by removing portions of the band edges in which the bandpass filter rolled off significantly, as well as channels that were identified as containing significant (remaining) RFI.The presence of RFI was determined by manually inspecting the average profile for each receiver as a function of frequency and pulse phase.
Average profiles in the sixty-day periods immediately before and after the event, as well as in the most recent available data, are shown in Figure 1.We compared these profiles with templates for each band, created by averaging the profiles from all available observations made before the event.Then, for each profile shown in Figure 1, we subtracted a best-fit scaled, aligned copy of the corresponding template.In the presence of a shape change, the correct alignment between the profile and the template is ambiguous, because the normal procedure for fitting TOAs breaks down, as discussed in more detail below.To illustrate this, in quantifying the degree of the shape change, we used two different methods of alignment.
In the "fixed-phase" method, the template is aligned by extrapolating the phase predicted by a timing model fit to data before the event.As a result, the reference point for comparison in this case is the best prediction based on data available before the change.In effect, this method assumes that the "true" TOA (based on the rotation of the pulsar) has not changed as a result of the event.For this assumption to be a reasonable one, changes in the rotational phase of the pulsar due to intrinsic spin noise have to be small enough to be neglected.This is much more likely to be the case for an MSP such as J1713+0747 than for a canonical pulsar.In the "fit-phase" method, on the other hand, the template was aligned by fitting for a TOA in the ordinary manner (i.e., using Fourier-domain matched filtering, as described by Taylor 1992), and shifting the template accordingly.This necessarily results in a smaller shape difference -in fact, the smallest possible shape difference for any choice of alignment -but there is no good reason to believe that the corresponding TOA is connected to the rotation of the pulsar.However, TOAs calculated in this way can serve as a useful point of comparison (see Section 2.4 below).
To quantify the degree to which the shape changed in each band, a fractional shape difference was derived for each combination of profile and alignment method.The results are shown in Table 2 below.The shape difference was derived by first aligning the template with the profile, using either the "fixed-phase" or "fit-phase" method, and then fitting for the scale factor such that the aligned and scaled template (i.e., the profile model) most closely matched the observed profile.The model was then subtracted from the profile to produce a profile residual.The shape difference was then calculated by dividing the root-mean-square (RMS) of the profile residual by the RMS of the profile model (calculated after normalizing the profile to unit amplitude).
The frequency dependence of the profile shape, both before and after the event, can be seen more clearly in Figure 2, which shows the frequency-dependent profiles across all three bands in the same sixty-day periods used in Figure 1. Figure 3 includes frequency-averaged profiles for the full set of GBT and CHIME data observed between February 21, 2020 and February 1, 2023, illustrating the time dependence of the profile shape.Because the frequency dependence of the profile shape is greatest in the CHIME band, Figure 4 is also included, demonstrating the time dependence of the profile shape in each of four 100 MHz sub-bands of the CHIME data.To quantify the degree to which the shape had changed, we produced profile residuals by subtracting a best-fit scaled copy of the template from each of the profiles in Figure 3.The results are shown in Figure 5.In producing the profile residuals shown in Figure 5, we used the "fixed-phase" method; i.e., the phase of the template subtracted from each profile was based on the timing model rather than fit to the data.As a result, the changes seen in Figure 5 are relative to expectations based on data from before the shape change.

Principal component analysis
We then performed principal component analysis (PCA) on the profile residuals.The vectors v ik in this basis are called the principal components.They are eigenvectors of the sample covariance matrix, C ij , which has the eigendecomposition Here n is the number of observations, and the indices i and j range over the m phase bins.There are r principal components, where r = min(m, n) is the rank of the covariance matrix, C ij .The principal components are orthonormal in the sense that (2) The eigenvalue, λ ℓ , corresponding to a principal component represents the amount of variance it describes in the data.Principal components are conventionally sorted in order of decreasing eigenvalue, and typically only the first few, most significant components are of interest.
Observations can be decomposed into principal components as The coefficients, z kℓ , in this expansion will be referred to as the principal component amplitudes.They represent the degree to which a particular principal component contributes to a particular observation, and their scale is determined by the normalization of the principal components.They can be computed from the data as which follows from the orthonormality of the principal components (equation 2).We performed PCA separately on the frequency-averaged profile residuals in each band, using a singular value decomposition of the profile residuals.The most significant principal components in each band (two in the GBT 1500 MHz band, and three in each of the GBT 820 MHz and CHIME bands) are shown in Figure 6, along with the correspoding amplitudes, z ij .No other principal components were significant.Multiplying the amplitudes by the corresponding principal components and adding them up, as in equation ( 3), produces a time-dependent model for the shape of the profile, which will be discussed further in Section 3 below.Uncertainties on the principal component amplitudes, z ij , were calculated by assuming the observed profile residuals, x ik , were subject to additive white noise.The amplitude of the noise in each profile was determined empirically from its discrete Fourier transform (DFT), by taking the standard deviation of the highest-frequency half of the Fourier components.Since the signal of interest is low-frequency, these highest frequency components consist almost entirely of noise, and, because the (appropriately normalized) DFT is a unitary transformation, the amplitude of this noise is directly related to the amplitude of the additive white noise in the time domain.The resulting uncertainty estimate for x ij was then propagated through equation 4, assuming z ij is held fixed, to calculate the uncertainty on z ij .This is a reasonably good estimate of the uncertainty in estimating the shape of the observed profile.However, it does not take into account any intrinsic variability in the true shape, which, as seen in section 3 below, may be significant.

TOA and DM estimation
Because standard TOA and DM estimates rely on the assumption that pulse shapes are consistent over time, they are necessarily biased in the presence of a pulse shape change.However, because such TOA and DM estimates are routinely produced as part of pulsar timing analyses, it is common for events like this one to be analyzed and quantified initially in terms of their effects on apparent TOA and DM.To facilitate comparison with such previous analyses, we produced TOA and DM estimates in the ordinary way; i.e., using the "fit-phase" procedure described above.As in the analyses described previously, the templates used were derived only from data collected prior to the onset of the shape change.This ensures that the TOA and DM estimates thus obtained, although unrelated to the true rotational phase of the pulsar or line-of-sight electron density, capture changes relative to the earlier, steady state.
More sophisticated TOA estimation techniques, which seek to correct for the effect of the shape change and recover the underlying rotation phase of the pulsar, are possible.For example, as shown by Demorest (2007) and further explored by Osłowski et al. (2011) in the context of pulse jitter, measured principal component amplitudes can be exploited to partially remove TOA estimation bias.Similar corrections can also be performed using other profile representations, including shapelet bases (Lentati et al. 2015;Lentati et al. 2017) and bases derived from single-pulse clustering algorithms (Kerr 2015).In this context, a caveat applies, in that if there is a true change in TOA associated with the pulse shape change and varying in amplitude together with it, such methods will also remove the TOA change.We intend to explore the use of these more sophisticated estimation techniques in future work, but here we consider only the naïve TOA and DM estimates described above.To produce these estimates, we used the software package PulsePortraiture (Pennucci et al. 2014;Pennucci 2019).We derived a frequency-dependent profile model in each band by aligning and averaging the data observed prior to the event, performing PCA on the resulting per-channel profiles, and fitting a B-spline model to the corresponding principal component amplitudes, using the method described by Pennucci et al. (2014).In the frequency band between 720 and 800 MHz, which is visible both to the GBT 820 MHz receiver and to CHIME, the two profile models are largely consistent, with a shape difference (in the sense of Table 2) of less than 1.1%.Using these models, we performed two types of fits to each observation.First, the full frequency-dependent profile model was fit to the data, producing a DM estimate and a single TOA estimate, referenced to the center frequency of the data.We refer to the TOA and DM estimates produced this way as "wideband" TOAs and DMs.Then, we divided each band into 12.5 MHz channels and fit profile models derived from the full frequency-dependent model to the data in each channel separately, producing estimated TOAs in each channel.We refer to the results as "narrowband" TOAs.TOA residuals were produced from the narrowband TOAs by subtracting the average TOA in each frequency channel prior to the event from all of the TOAs in that channel, to eliminate the effect of any frequency dependence present prior to the shape change.Additionally, wideband DM residuals were produced by subtracting the mean DM in each receiver band prior to the shape change from the wideband DM estimates in that band.The TOA and DM residuals in each receiver band are shown in Figure 7, and the TOA residuals are shown broken down by frequency in Figure 8.As mentioned above, the resulting TOA and DM measurements are highly covariant with the pulse shape, and should be understood as demonstrating how the event would appear in a conventional pulsar timing analysis, and not as describing true accompanying changes in the pulsar's spin or the column density of electrons along the line of sight.

RESULTS
The average profiles seen in Figure 1 illustrate the form of the shape change at different radio frequencies.As indicated in Table 2, all bands are affected to a similar degree: the shape difference is approximately 40 to 60 percent, compared to a fixed-phase template, or 20 to 40 percent compared to a template fit to the data in phase.Regardless of the alignment method used, this is much larger than the profile shape changes reported to accompany the 2016 event in this pulsar by Lin et al. (2021), which were approximately 1 to 4 percent of the peak profile amplitude.
Because the profile residuals were calculated relative to a fixed-phase template, any phase shift associated with the shape change event would have been absorbed into them, in the form of a shape change proportional to the derivative of the template shape (cf.Osłowski et al. 2011, section 4).Other kinds of shape changes may coincidentally align partially with the derivative of the template shape, but will in general have no relation to it.For the GBT 1500 MHz, GBT 820 MHz, and CHIME bands respectively, we find that 79.4%, 85.9%, and 36.6% of the variance of the first principal component is in the direction of the derivative of the template shape, meaning that no more than roughly this fraction of the variance could be due to an associated phase shift.
Both before and after the onset of the event, the profile shows a noticeable dependence on frequency, as can be seen in Figure 2. Five main emission components can be identified: in addition to the bright central peak (C), there are two .Average profiles before, during and after the event, as observed at GBT and CHIME.The left and center profiles are averaged over sixty-day periods immediately before (MJDs 59261-59321) and after the shape change (MJDs 59321-59381), while the right profile is averaged over a sixty-day period 595 days later (MJDs 59916-59976), after the pulse shape had mostly recovered.Each sixty-day period includes two GBT observations in each band and almost 60 CHIME observations.Profiles are averaged across the bands described in Table 1.For the GBT data, the total linear polarization (L) and circular polarization (V ) are shown in addition to the total intensity (I), while, for the CHIME data, which are not polarization calibrated, only the total intensity is shown.
leading components (A, B) and two trailing components (D, E).One of these (D) is evident in the original profile only as a slight change in slope on the trailing edge of the main peak, together with a feature in the polarized emission, but appears much more distinct after the shape change.After the shape change, the main peak (C) is narrower, while the leading "shoulder" component (B) is brighter, and the trailing shoulder (D) is more pronounced.The second trailing component (E) also appears to be brighter than previously, when compared to the other components.Notably, the overall amplitude of emission from the pulsar is expected to vary between observing epochs as a result of scintillation, so it is much more difficult to measure absolute changes in the brightness of particular components than it is to measure relative changes in amplitude between the components.Figure 3 demonstrates the fact that the profile has recovered toward its original shape over the months that have elapsed since the event.At 1.5 GHz, the leading shoulder component is brightest, and the trailing shoulder most distinct, immediately after the shape change, with both effects tending to diminish over time.At 820 MHz, the   leading shoulder appears to have briefly become approximately as bright as the main peak, and gradually declines in brightness over time, while, as at 1.5 GHz, the trailing shoulder blends more and more with the main pulse.The recovery of the profile toward its original state is also clear in Figure 5, which shows the profile residuals.The residuals can be seen to gradually decrease over time, while approximately maintaining their shape.As seen in Figure 6, this recovery can be described primarily as a change in the amplitude of the first principal component, z 1 , as a function of  1, and normalized to a constant peak intensity.Both GBT and CHIME observations are shown as horizontal bands whose lower edge corresponds to the observation date.For the CHIME observations, the band height is always one day, so gaps appear on days when CHIME did not observe the pulsar.For the much less frequent GBT observations, the bands are extended vertically until the date of the next observation, or, in the case of the last observation, for an additional 30 days.This means that the shape change appears to occur slightly earlier at CHIME than at GBT, when in fact, the lower observing cadence at GBT makes this impossible to determine.
epoch.The form and timescale of the recovery can be quantified by fitting a parametric model to z 1 (t).We fit three such models to the data in each band: a decaying exponential (Model E): a decaying exponential with persistent offset (Model O): and a power-law model (Model P): The results can be seen in Table 3.In each case, t 0 was taken to be MJD 59321.0exactly; the true onset of the shape is uncertain, but differs from this by no more than about 12 hours.Data from before this point were excluded from the fit.Model E, the simplest, did not produce a satisfactory fit, but models O and P fit approximately equally well.Even for these better-fitting models, the resulting reduced chi-squared statistic was very large, reflecting the fact that there is additional short-time variability not accounted for by the nominal uncertainties.This is perhaps not surprising, since these uncertainties include only estimation error arising from white noise in the profile residuals, and not any intrinsic variability in pulse shape.Nevertheless, models O and P appear to describe the long-term evolution of the first principal component amplitude in each band relatively well.Note that model O predicts that some fraction of the shape change will remain indefinitely, whereas model P predicts that the profile shape will continue to recover, albeit more slowly.Unfortunately, since the data do not exhibit a clear preference between the two models, it is not possible to draw firm conclusions as to which of these predictions is more likely to hold.
The effect of the shape change on TOA estimates depends on frequency in a non-monotonic fashion, as demonstrated in Figure 7.In the GBT 820 MHz and 1.5 GHz data, the effect of the shape change is consistently larger at lower frequencies.This frequency dependence manifests itself as an apparent change in DM peaking at approximately −7.5 × 10 −3 pc cm −3 in the 1.5 GHz band and −2.5 × 10 −3 pc cm −3 in the 820 MHz band, an order of magnitude larger than that seen in the previous two, smaller chromatic timing events, but in the same direction.This behaviour is consistent with what would be expected if the event were caused by the passage of an under-dense region of the ISM through the line of sight, which is the explanation preferred by Lam et al. (2018) and Lin et al. (2021) for the two previous chromatic timing events.However, the CHIME data indicate that, at the lowest frequencies, this frequency dependence reverses itself.This means that the apparent change in DM measured in the CHIME band has the opposite sense, peaking at approximately 2.5 × 10 −3 pc cm −3 , something which is difficult to explain in a model which relies on changes to the ISM along the line of sight to the pulsar.
Importantly, conventional TOA and DM measurements like those made here are always covariant with pulse shape changes, which makes it impossible to measure an absolute TOA without making assumptions about the pulse shape, and similarly impossible to measure an absolute DM without making assumptions about the frequency dependence of the pulse shape.Because we make the deliberately naïve choice to use a profile model based only on data prior to the shape change, the TOA and DM measurements given in this section and in Figures 7 and 8 should be understood as an empirical means of characterizing the shape change, rather than as a reflection of real changes to the underlying "true" TOA (defined with respect to the rotation of the pulsar) or DM (i.e., the column density of electrons along the line of sight)., regions where the observed intensity is greater than that predicted by the model, and blue represents negative residuals, regions where the observed intensity is less than the model prediction.Also, as in Figure 3, the shape change appears to occur slightly earlier at CHIME than at GBT, but this is merely the result of the lower observing cadence at GBT.
In many respects, the current event resembles a larger version of the two previous chromatic timing events.All three correspond to abrupt decreases in apparent DM, at least at frequencies above approximately 800 MHz, and all three show at least some evidence for associated profile shape changes, which are only gradually frequency-dependent.The decay time of the new event is approximately 156 days at L-band (see Table 3), which can be compared with the 62 and 25 day decay times derived for the previous two events by Lam et al. (2018).While the new event takes significantly longer to decay away, the decay time is of a similar order of magnitude.It may therefore be reasonable to assume that the new event has the same physical origin as those previous events.
The shape change is too strongly correlated in frequency and time to be explained by diffractive interstellar scintillation (DISS).DISS can lead to apparent changes in pulse shape when observations of a pulsar with a frequencydependent pulse shape are integrated across a wide band.This happens because scintillation makes the signal at certain frequencies appear brighter than at other frequencies, which biases the shape of the integrated pulse toward the shape at the brightened frequencies.However, this can be ruled out as the cause of the shape changes seen here.For J1713+0747, the bandwidth, ∆ν d , and time, ∆t d , over which the scintillation pattern decorrelates, have been measured to be approximately 22 MHz and 48 min, respectively (Levin et al. 2016;Turner et al. 2021).As described in Section 3 above, the observed shape change is correlated over much larger bandwidths (hundreds of MHz) and times (months) than this, and its magnitude exceeds the degree of profile variation with frequency observed prior to the event.Moreover, the shape change can be also be seen in subbands narrower than the 22 MHz scintillation bandwidth, in a form similar to what is seen in the integrated profiles. .Principal component analysis of the profile residuals in each band.For each of the three bands described in Table 1 (GBT 1500 MHz, GBT 820 MHz, CHIME 600 MHz), the three most significant principal components are shown at left, along with the average total intensity profile in that band.To the right, the corresponding principal component amplitudes are shown as a function of time.Adding the principal components, weighted according to their respective amplitudes, to the average profile gives a model of the profile shape at any particular observing epoch.  .Apparent band-averaged TOA (top) and DM (bottom) residuals before and after in the period around the shape change event, as measured using a profile model fit to data from before the event.Both the TOA and DM measurements are highly covariant with the pulse shape, and so these should be understood as demonstrating how the event would appear in a conventional pulsar timing analysis, and not as describing true accompanying changes in the pulsar's spin or the column density of electrons along the line of sight.The fact that the TOA residuals are a non-monotonic function of frequency, reaching their most negative value at around 750 MHz (cf. Figure 8), can be seen in both of these plots: the amplitude of the TOA residuals is greater in the 820 MHz band than at either 600 MHz or 1500 MHz; and the apparent DM residuals, which roughly represent the slope of the TOA residuals as a function of frequency, are positive at 600 MHz and negative at 1500 MHz.
An event possibly analogous to the profile shape change described here was observed in the MSP J1643−1224 beginning around MJD 57080 (February 27, 2015;Shannon et al. 2016).Like the recent event in J1713+0747, the 2015 event in J1643−1224 has a rise time of no more than a few days, and subsequently decays over the course of several months.Shannon et al. (2016) attribute this event to a magnetospheric disturbance of unspecified nature.It also shows a significant dependence on radio frequency, but one which is inverted compared to that seen in the first two events in J1713+0747, with more pronounced shape changes, and correspondingly larger timing effects, seen at higher frequencies.In observations with the 64-meter Parkes radio telescope (Murriyang), Shannon et al. (2016) found that the change in the shape of J1643−1224's profile was most prominent at 10 cm (3 GHz), less noticeable at 20 cm (1.5 GHz), and undetectable at 50 cm (600 MHz).Notably, however, Brook et al. (2018) examined GBT observations of J1643−1224 taken during the same time period, and found that the effect was stronger at 800 MHz than at 1.5 GHz.This may point to a non-monotonic frequency dependence for the J1643−1224 event, as is seen here for the recent J1713+0747 event.
Several different hypotheses may be entertained as to the origin of the shape change event.One possibility is that the event might be due to lensing of the radio emission by a discrete structure in the ionized ISM along the line of sight.This is the generally accepted explanation for the 1997 event observed in the Crab pulsar (Backer et al. 2000;Graham Smith et al. 2011), and is the explanation proposed by Lam et al. (2018) and Lin et al. (2021) for the previous two events seen in PSR J1713+0747.Plasma lensing is expected to produce radio frequency-dependent changes in TOA measurements, and can also produce multiple images which are delayed in time by amounts comparable to the width of the pulse (0.1-1 ms), creating apparent profile shape changes (Lam et al. 2016;Cordes et al. 2017).However, the complex morphology of the shape change, which includes a sharpening of parts of the trailing edge, is not easily explained as a combination of a small number of lensed images.The non-monotonic frequency dependence of the TOA residuals may also point away from an ISM interpretation, but it is far from a perfect indicator.Simple changes in DM produce timing delays that are proportional to ν −2 , and, while scattering and frequency-dependent DM effects can modify this dependence, they typically do so only modestly.However, ray tracing through a plasma lens can produce more complex frequency-dependent effects: refraction occurs due to the gradient of DM across the lens, and caustics can be produced.Cordes et al. (2017) studied these phenomena in the context of fast radio bursts, finding that a variety of effects were possible, including cases where the frequnency dependence of arrival times differs markedly from the ν −2 scaling expected from dispersion alone.While plasma lensing does not seem to be the most likely explanation for the present event, it cannot be ruled out entirely without a more detailed analysis.Some previously observed profile shape changes have been caused by geodetic precession, but this can be ruled out for J1713+0747.In binary neutron stars such as PSR B1913+16 and the double pulsar, PSR J0737−3039A/B, geodetic precession can cause changes in the viewing geometry, leading to pulse shape changes (Weisberg et al. 1989;Kramer 1998;Breton et al. 2008).However, J1713+0747 has a white dwarf companion in a wide orbit with a 67.8day period (Splaver et al. 2005), which places it too far away for geodetic precession to be significant (the expected precession rate is a mere 2.2 mas/yr).Additionally, for geodetic precession to produce detectable effects, a pulsar's spin direction must be significantly misaligned with its orbit, and the process of mass transfer that leads to the formation of MSP-white dwarf binaries like the J1713+0747 system tends to drive spin-orbit alignment.In any case, the rapid onset of the shape change seen here is inconsistent with such an origin.
Many other pulsars exhibit mode-changing, in which the profile abruptly changes back and forth between 2 or more shapes (e.g.Backer 1970, Wang et al. 2021, Miles et al. 2022).This phenomenon has been shown to extend to millisecond pulsars (Mahajan et al. 2018).Gradual recovery of the profile toward a previous state, as is seen in this case, has not previously been observed in mode-changing pulsars.However, in some cases, there is a pronounced relationship between the profile state and the spin evolution of the pulsar (Kramer et al. 2006;Lyne et al. 2010;Shaw et al. 2022).In the case of PSR B1828−11, the fraction of time spent in each of the two profile modes varies roughly periodically with time and is directly related to the local value of the spin-down rate (Stairs et al. 2019).Some of these "magnetospheric switching" pulsars can have very long-lived runs in a given profile state (e.g.PSR B2035+36, Shaw et al. 2022).It is conceivable that the profile change observed here in PSR J1713+0747 is seen to "decay" back to its pre-event state because of slow changes in the fraction of time spent in each of two extreme modes.Detailed inspection of high-time-resolution data will be needed to determine whether or not this is the case.
In some cases, profile shape changes might be expected to accompany pulsar glitches.Glitches are a type of abrupt spin-up event, thought to be symptoms of the rapid transfer of angular momentum from a pulsar's interior to its crust through the unpinning of superfluid vortices (Anderson & Itoh 1975;Haskell & Melatos 2015).Although they are most often seen in young pulsars, two glitches have been observed in MSPs (Cognard & Backer 2004;McKee et al. 2016).Glitches are sometimes followed by a period of months in which the pulsar's period and period derivative recover toward their pre-glitch values.The rapid onset and subsequent recovery over seen in the shape change event in J1713+0747 are somewhat reminiscent of this phenomenon.However, glitches are not usually assoicated with profile shape changes, and in the few cases where profile changes associated with glitches have been observed, this is not what they look like.Weltevrede et al. (2011) observed additional intermittent pulse components in the aftermath of a glitch in PSR J1119−6127, and Keith et al. (2013) observed a relationship between glitching and emission state changes in PSR J0742−2822.More recently, Palfreyman et al. (2018) managed to observe the Vela pulsar, PSR B0833−45, continuously over a period including the onset of a glitch, and saw changes in a few single pulses at the moment the glitch occurred.None of these scenarios is a particularly good analog for what is seen here in J1713+0747.The primary feature of a glitch is an achromatic step in the pulsar's spin frequency, which is not present in the current case.Furthermore, the shape changes which have previously been seen to accompany glitches are either transient or intermittent, rather than being sustained over several months with gradual decay, as seen here.
A more exotic possibility is that the shape change resulted from the injection of an asteroid into the pulsar's magnetosphere.Brook et al. (2014) argue that a new emission component seen in the canonical pulsar PSR B0736−40 in 2006 (Karastergiou et al. 2011) was caused by such an asteroid incursion.Asteroidal material which entered the pulsar's magnetosphere, altering the flow of particles and thus the radio emission, could result in observable changes to a pulsar's radio emission (Cordes & Shannon 2008).However, this is relatively unlikely for MSPs such as PSR J1713+0747.Since the light-cylinder radius r LC = cP is 100 to 1000 times smaller for MSPs than for CPs, an asteroid is much more likely to be destroyed before entering the magnetosphere (Cordes & Shannon 2008).
Overall, while the origin of the shape change event in J1713+0747 remains far from clear, it seems most likely that it was caused by some kind of change that took place within the pulsar's magnetosphere, as ISM-based explanations cannot easily account for its complex morphology.It is sufficiently different from a glitch that it must be considered a separate phenomenon.If it indeed has the same origin as the two previous events seen in J1713+0747 and the 2015 event in J1643−1224, it is by far the largest such event ever observed, and may reveal important aspects of the general nature of events of this kind.

CONCLUSIONS
The recent shape change observed in PSR J1713+0747 most likely originated in the pulsar's magnetosphere.In several ways, it resembles a larger version of the two previously observed chromatic timing events (Demorest 2007;Lam et al. 2018;Lin et al. 2021), which have previously been attributed to lensing of the pulsar emission by underdense regions in the ISM.As in these past events, the DM, as measured between 820 and 1500 MHz, has decreased abruptly, and appears to be recovering toward its original value on a timescale of several months.Although it is frequencydependent, the event shows non-monotonic behavior inconsistent with a simple change in DM, and is accompanied by pulse shape changes that are nearly achromatic.There is some evidence that the same may be true of the two previous events.
ISM propagation effects in pulsar timing typically produce chromatic changes in TOAs.However, it is difficult to produce complex shape changes in this way.Lensing may produce multiple images of the pulse that interfere with each other to produce an altered profile shape, but systematic changes in the widths and relative amplitudes of pulse components, such as those seen in the event described here, would be much more likely given a magnetospheric origin.
The new event bears some resemblance to the profile shape change seen in PSR J1643−1224 in early 2015 (Shannon et al. 2016), which is thought to have had a magnetospheric origin.In that case, the frequency dependence was inverted, with the shape changes being stronger at higher frequencies.The non-monotonic frequency dependence of the event described here, however, means that an analogous origin for the J1643 event is very likely.An event of similar form may also have occurred in PSR J1640+2224 in mid-2012 (Hazboun et al., in prep.).Events like these may represent an entirely new phenomenon that will complicate millisecond pulsar timing, including pulsar timing array searches for gravitational waves.Fully accounting for them will likely require the use of time-varying templates referenced to the same fiducial pulse phase, as was done for PSR J0737−3039B, the young pulsar component of the double pulsar, whose pulse changes shape as a result of geodetic precession, in Noutsos et al. (2020).We intend to explore the use of time-varying templates in more detail in future work.
The trend toward increasing fractional bandwidths in pulsar observing has created a need for pulsar timing techniques which can account for profile shape variations across a wide band, without imposing an undue computational burden on later analysis.Compared to splitting a wide band into many narrow subbands and generating traditional narrowband TOAs within each subband, wideband techniques like those of Pennucci (2019) can account for such shape variations in a natural way while reducing the amount of data that must be processed in subsequent analysis steps to a single TOA and DM (rather than a TOA in each subband) per epoch.The downside of this approach is that other kinds of frequency dependence, such as time-variable scattering or the kind of time-and-frequency dependent event seen here, can no longer be dealt with at the level of TOA residuals, since information about frequency dependence within each band, beyond the DM, is lost.There is a potential solution to this problem, however, which is to deal with these frequency-dependent effects directly at the level of the profiles, by adding appropriate terms to the wideband TOA log-likelihood (cf.Pennucci 2019;Alam et al. 2021b).Indeed, such an approach may be the only way to deal with an event like the one described here within the a fully wideband framework.
Events such as the present one will present a challenge for low-frequency gravitational wave searches by PTAs.However, there are good reasons to believe this challenge is surmountable.Gravitational wave signals are expected to be dominated by very low-frequency components, appear independent of radio frequency, and show a characteristic quadrupolar pattern of spatial correlations between pulsars, quantified by the Hellings & Downs curve (Hellings & Downs 1983;Cornish & Sesana 2013;Arzoumanian et al. 2020).Because of its limited extent in time (in particular, data taken before the onset of the event is unaffected by it), its dependence on radio frequency, and its appearance in only a single pulsar, this event is not likely to be strongly covariant with gravitational-wave signals, and it should in principle be possible to model and subtract its effects without severely affecting the gravitational-wave sensitivity of a PTA data set.
Figure1.Average profiles before, during and after the event, as observed at GBT and CHIME.The left and center profiles are averaged over sixty-day periods immediately before (MJDs 59261-59321) and after the shape change (MJDs 59321-59381), while the right profile is averaged over a sixty-day period 595 days later (MJDs 59916-59976), after the pulse shape had mostly recovered.Each sixty-day period includes two GBT observations in each band and almost 60 CHIME observations.Profiles are averaged across the bands described in Table1.For the GBT data, the total linear polarization (L) and circular polarization (V ) are shown in addition to the total intensity (I), while, for the CHIME data, which are not polarization calibrated, only the total intensity is shown.

Figure 2 .
Figure2.Pulse profiles as a function of frequency and phase before, during, and after the shape change, averaged over the same sixty-day periods as in Figure1.All three portraits are de-dispersed with the same dispersion measure, 15.9638 pc cm −3 .Observations from all three receivers are included, with CHIME data plotted over GBT 820 MHz data where the bands overlap.A distinct leading component appearing only at the lowest frequencies is evident.

Figure 3 .
Figure 3. of the shape change event as seen by the GBT and CHIME.Profiles observed at GBT, using the VEGAS backend, are seen in the left two panels, and those observed at CHIME are shown in the right panel.All profiles are averaged across the bands described in Table1, and normalized to a constant peak intensity.Both GBT and CHIME observations are shown as horizontal bands whose lower edge corresponds to the observation date.For the CHIME observations, the band height is always one day, so gaps appear on days when CHIME did not observe the pulsar.For the much less frequent GBT observations, the bands are extended vertically until the date of the next observation, or, in the case of the last observation, for an additional 30 days.This means that the shape change appears to occur slightly earlier at CHIME than at GBT, when in fact, the lower observing cadence at GBT makes this impossible to determine.

Figure 4 .
Figure 4.The shape change event as seen in four 100 MHz sub-bands of the CHIME data.The low-frequency component seen in Figure 2 is visible in the lower three sub-bands, and the high frequency component is visible in the upper three.

Figure 5 .
Figure 5. Profile residuals around the shape change event, derived from the profiles in Figure 3 by subtracting a "fixedphase" template (cf.Section 2.2).As in Figure 3, GBT (VEGAS) observations are shown in the left two panels, and CHIME observations are shown in the right panel.Red represents positive residuals; i.e., regions where the observed intensity is greater than that predicted by the model, and blue represents negative residuals, regions where the observed intensity is less than the model prediction.Also, as in Figure3, the shape change appears to occur slightly earlier at CHIME than at GBT, but this is merely the result of the lower observing cadence at GBT.
Figure6.Principal component analysis of the profile residuals in each band.For each of the three bands described in Table1(GBT 1500 MHz, GBT 820 MHz, CHIME 600 MHz), the three most significant principal components are shown at left, along with the average total intensity profile in that band.To the right, the corresponding principal component amplitudes are shown as a function of time.Adding the principal components, weighted according to their respective amplitudes, to the average profile gives a model of the profile shape at any particular observing epoch.

Figure 7
Figure7.Apparent band-averaged TOA (top) and DM (bottom) residuals before and after in the period around the shape change event, as measured using a profile model fit to data from before the event.Both the TOA and DM measurements are highly covariant with the pulse shape, and so these should be understood as demonstrating how the event would appear in a conventional pulsar timing analysis, and not as describing true accompanying changes in the pulsar's spin or the column density of electrons along the line of sight.The fact that the TOA residuals are a non-monotonic function of frequency, reaching their most negative value at around 750 MHz (cf.Figure8), can be seen in both of these plots: the amplitude of the TOA residuals is greater in the 820 MHz band than at either 600 MHz or 1500 MHz; and the apparent DM residuals, which roughly represent the slope of the TOA residuals as a function of frequency, are positive at 600 MHz and negative at 1500 MHz.

Figure 8 .
Figure 8. Dependence of the apparent TOA residuals on frequency.The left panel shows the average TOA residual as a function of frequency, in each of two thirty day periods: immediately after the shape change (MJD 59321-59351, marked by dots), and at the end of the data presented here (MJD 59621-59651, marked by crosses).The right panel shows the TOA residual as a function of both time and frequency, averaged into 25 MHz × 30 day bins.The frequency dependence is not consistent with the ν −2 form expected for dispersive behavior, but instead is non-monotonic, with the TOA residuals reaching their most negative value at around 750 MHz.

Table 1 .
Summary of observations The average cadence of GBT 1500 MHz observations was reduced to equal that of the 820 MHz observations after March 2021 due to the end of the NANOGrav high-cadence observing program.

Table 2 .
Degree of shape change in each band.Fractional shape differences are calculated using the root-mean-square method described in Section 2.2.

Table 3 .
Time evolution of first principal component