An unusual pulse shape change event in PSR J1713+0747 observed with the Green Bank Telescope and CHIME
0000-0003-1082-2342]Ross J. Jennings
NANOGrav Physics Frontiers Center Postdoctoral Fellow
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA
0000-0002-4049-1882]James M. Cordes
0000-0002-2878-1502]Shami Chatterjee
Department of Astronomy, Cornell University, Ithaca, NY 14853, USA
Cornell Center for Astrophysics and Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY 14853, USA
0000-0001-7697-7422]Maura A. McLaughlin
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA
0000-0002-6664-965X]Paul B. Demorest
National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA
Zaven Arzoumanian
X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA
0000-0003-2745-753X]Paul T. Baker
Department of Physics and Astronomy, Widener University, One University Place, Chester, PA 19013, USA
0000-0003-4046-884X]Harsha Blumer
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA
0000-0003-3053-6538]Paul R. Brook
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA
0000-0001-7587-5483]Tyler Cohen
Deptartment of Physics, New Mexico Institute of Mining and Technology, 801 Leroy Place, Socorro, NM 87801, USA
0000-0002-2578-0360]Fronefield Crawford
Department of Physics and Astronomy, Franklin & Marshall College, P.O. Box 3003, Lancaster, PA 17604, USA
0000-0002-6039-692X]H. Thankful Cromartie
NASA Hubble Fellowship: Einstein Postdoctoral Fellow
Cornell Center for Astrophysics and Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY 14853, USA
0000-0002-2185-1790]Megan E. DeCesar
George Mason University, resident at the Naval Research Laboratory, Washington, DC 20375, USA
0000-0001-8885-6388]Timothy Dolch
Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA
Eureka Scientific, Inc., 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA
0000-0001-7828-7708]Elizabeth C. Ferrara
Department of Astronomy, University of Maryland, College Park, MD 20742, USA
Center for Research and Exploration in Space Science and Technology, NASA/GSFC, Greenbelt, MD 20771, USA
Goddard Space Flight Center, Greenbelt, MD 20771, USA
0000-0001-8384-5049]Emmanuel Fonseca
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA
0000-0003-1884-348X]Deborah C. Good
Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046, USA
Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
0000-0003-2742-3321]Jeffrey S. Hazboun
University of Washington Bothell, 18115 Campus Way NE, Bothell, WA 98011, USA
0000-0001-6607-3710]Megan L. Jones
Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA
0000-0001-6295-2881]David L. Kaplan
Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA
0000-0003-0721-651X]Michael T. Lam
School of Physics and Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA
Laboratory for Multiwavelength Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA
T. Joseph W. Lazio
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA
0000-0003-1301-966X]Duncan R. Lorimer
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA
0000-0001-5373-5914]Jing Luo
Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON M5S 3H4, Canada
0000-0001-5229-7430]Ryan S. Lynch
Green Bank Observatory, P.O. Box 2, Green Bank, WV 24944, USA
0000-0002-2885-8485]James W. McKee
Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 Saint George Street, Toronto, ON M5S 3H8, Canada
E.A. Milne Centre for Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK
Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK
0000-0003-2285-0404]Dustin R. Madison
Department of Physics, University of the Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA
0000-0001-8845-1225]Bradley W. Meyers
Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
0000-0002-4307-1322]Chiara M. F. Mingarelli
Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046, USA
0000-0002-6709-2566]David J. Nice
Department of Physics, Lafayette College, Easton, PA 18042, USA
0000-0001-5465-2889]Timothy T. Pennucci
Institute of Physics, Eötvös Loránd University, Pázmány P.s. 1/A, 1117 Budapest, Hungary
0000-0002-8509-5947]Benetge B. P. Perera
Arecibo Observatory, University of Central Florida, HC3 Box 53995, Arecibo, PR 00612, USA
0000-0002-8826-1285]Nihan S. Pol
Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235, USA
0000-0001-5799-9714]Scott M. Ransom
National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA
0000-0002-5297-5278]Paul S. Ray
U.S. Naval Research Laboratory, Washington, DC 20375, USA
0000-0002-7283-1124]Brent J. Shapiro-Albert
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA
Giant Army, LLC, 915A 17th Ave, Seattle, WA 98122, USA
0000-0002-7778-2990]Xavier Siemens
Department of Physics, Oregon State University, Corvallis, OR 97331, USA
Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA
0000-0001-9784-8670]Ingrid H. Stairs
Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
0000-0002-1797-3277]Daniel R. Stinebring
Department of Physics and Astronomy, Oberlin College, Oberlin, OH 44074, USA
0000-0002-1075-3837]Joseph K. Swiggum
NANOGrav Physics Frontiers Center Postdoctoral Fellow
Department of Physics, Lafayette College, Easton, PA 18042, USA
0000-0001-7509-0117]Chia Min Tan
Department of Physics, McGill University, 3600 rue University, Montréal, QC H3A 2T8, Canada
McGill Space Institute, McGill University, 3550 rue University, Montréal, QC H3A 2A7, Canada
0000-0003-0264-1453]Stephen R. Taylor
Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235, USA
0000-0003-4700-9072]Sarah J. Vigeland
Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA
0000-0002-6020-9274]Caitlin A. Witt
Center for Interdisciplinary Exploration and Research in Astrophysics, Northwestern University, 1800 Sherman Ave, Evanston, IL 60201
Adler Planetarium, 1300 S. DuSable Lake Shore Drive, Chicago, IL 60605

The millisecond pulsar J1713+0747 underwent a sudden and significant pulse shape change between April 16 and 17, 2021 (MJDs 59320 and 59321). Subsequently, the pulse shape gradually recovered over the course of several months. We report the results of continued multi-frequency radio observations of the pulsar made using the Canadian Hydrogen Intensity Mapping Experiment (CHIME) and the 100-meter Green Bank Telescope (GBT) in a three-year period encompassing the shape change event, between February 2020 and February 2023. As of February 2023, 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 TOA residuals display a strong non-monotonic dependence on radio frequency, demonstrating that the event is neither a glitch (the effects of which should be independent of radio frequency, $\nu$) nor a change in dispersion measure (DM) alone (which would produce a delay proportional to $\nu^{-2}$). However, it does bear some resemblance to the two previous “chromatic timing events” observed in J1713+0747 [Demorest et al., 2013, Lam et al., 2016], as well as to a similar event observed in PSR J1643$-$1224 in 2015 [Shannon et al., 2016].

§ 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., 2022) and the sharp features in its pulse profile (effective width 0.54ms at L-band; based on data from Alam et al., 2021) 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., 2021; Alam et al., 2021), 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-hundred-meter Aperture Spherical Telescope (FAST) in China [Nan & Li, 2013, Hobbs et al., 2019] and the MeerKAT telescope in South Africa [Spiewak et al., 2022]. 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.3ns 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 <cit.>.

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 <cit.>. 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.1GHz) and L band (1228–1740MHz), respectively, as well as with the Kunming 40-meter telescope at S band (2150–2400MHz) and with FAST at L band (1050–1450MHz). 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.3e-3pc 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 $\nu^{-2}$ scaling expected from dispersion.

Furthermore, 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 frequency-dependent 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, $\Delta\nu_d$, and time, $\Delta t_d$, over which the scintillation pattern decorrelates, have been measured to be approximately 22MHz and 48min, respectively [Levin et al., 2016, Turner et al., 2021]. As described in Section <ref> below, 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 22MHz scintillation bandwidth, in a form similar to what is seen in the integrated profiles.

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 $\SI{-6e-4}{pc\,cm^{-3}}$, while for the second event, it is about $\SI{-4.3e-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.

Some possibly analogous events have been seen in other pulsars as well. Of particular relevance is the profile shape change 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. 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. Shannon et al., 2016 attribute the event to a magnetospheric disturbance of unspecified nature. While J1643$-$1224 is the only other MSP where similar pulse shape changes have been observed, a similar event, involving the appearance of a new emission component, was seen in the canonical pulsar B0736$-$40 in 2006 [Karastergiou et al., 2011]. Brook et al., 2014 attribute this event to an asteroid entering the pulsar's magnetosphere.

Another potential cause of pulse shape changes is 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.8-day 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.2mas/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.

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 <ref> 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 <ref>, 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 <ref>, we analyze the main features present in the data. Then, in Section <ref>, we consider several possible astrophysical interpretations for the observed phenomena. Finally, in Section <ref>, we summarize our results and conclude.

§ METHODS

§.§ Observation

l c c c c c c c

Summary of observations
Receiver     2cFrequency (MHz)     2cMJD     Number     Cadence (d$^{-1}$)     Median S/N

low     high     start     end
GBT 1500 MHz     $1100$     $1900$     58334     59945     136     0.084a     1410

GBT 820 MHz     $720$     $920$     58637     59916     37     0.029     938

CHIME 600 MHz     $400$     $800$     58596     59976     1096     0.794     65.6

aThe 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.

The observations considered here were made using the GBT and CHIME telescopes. The GBT observations were made using two receivers: the 800MHz band of the PF1 prime focus receiver (PF1/800), which has a bandwidth of 200MHz centered on 820MHz, and the L-band Gregorian focus receiver, which has 800MHz of bandwidth centered on 1500MHz. 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.5625MHz wide, with 512 such channels across the band at 1500MHz and 128 at 820MHz. GBT observations at both 820 and 1500MHz were made roughly monthly as part of the main NANOGrav observing program, and typically have integration times of about 30 minutes. Additional 1500MHz 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 1500MHz observations is lower after the onset of the shape change than before it.

The CHIME observations were made using the telescope's 256-antenna receiver system, which has 400MHz of bandwidth centered on 600MHz [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. 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.

A summary of all observations used in this paper is given in Table <ref>. 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 820MHz receiver and 70 observations with the 1500MHz receiver during this period. Of these, 16 of the 820MHz observations and 21 of the 1500MHz 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 prepared[The NANOGrav 15-year data set includes observations taken between July 30, 2004 (MJD 53216) and August 11, 2020 (MJD 59072) by two sets of pulsar backends: the Astronomical Signal Processor (ASP) at Arecibo and Green Bank Astronomical Signal Processor (GASP) at the GBT [Demorest], 2007], and the later Green Bank Ultimate Pulsar Processing Instrument (GUPPI) at the GBT [DuPlain et al., 2008] and its sibling instruments PUPPI at Arecibo and YUPPI at the Very Large Array. GUPPI was decommissioned in April 2020 and replaced by the VEGAS backend. Some VEGAS observations used in this paper were made simultaneously with GUPPI observations that are included in the 15-year data set, but the VEGAS observations are not included in the 15-year data set, nor are the GUPPI observations used here. PUPPI observations ceased after a cable failure at Arecibo in August 2020, and the Arecibo telescope ultimately collapsed on December 1, 2020, meaning that they were never resumed..  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[<https://doi.org/10.5281/zenodo.7236459>].

§.§ Post-processing and alignment

To prepare them for analysis, the GBT data were processed using NANOGrav's  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.9638pc cm^-3 and aligned using a pulsar ephemeris fit to the NANOGrav 15-year data. The exact parameter values used can be found in the Zenodo dataset. 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 <ref>. These average profiles were used to quantify the degree of the shape change event by fitting and subtracting templates based on the data before the event. 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. 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 <ref> 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 <ref> 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.

The frequency dependence of the profile shape, both before and after the event, can be seen more clearly in Figure <ref>, which shows the frequency-dependent profiles across all three bands in the same sixty-day periods used in Figure <ref>.  Figure <ref> 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 <ref> is also included, demonstrating the time dependence of the profile shape in each of four 100MHz 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 <ref>. The results are shown in Figure <ref>. In producing the profile residuals shown in Figure <ref>, 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 <ref> 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. PCA is a statistical method of dimensionality reduction, which linearly transforms a set of high-dimensional observations, $x_{ik}$ (here, the profile residuals), into the basis in which the greatest amount of variance can be described by the smallest number of variables.

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

\begin{equation}C_{ij} = \frac1{n}\sum_{k=1}^nx_{ik}\mkern2mu x_{jk} = \sum_{\ell=1}^r \lambda_\ell\mkern2mu v_{i\ell}\mkern2mu v_{j\ell}.
\end{equation}

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=\operatorname{min}(m, n)$ is the rank of the covariance matrix, $C_{ij}$. The principal components are orthonormal in the sense that

\begin{equation}\label{eqn:orthonormal-pcs}
\sum_{i=0}^m v_{i\ell}\mkern2mu v_{i\ell'} = \delta_{\ell\ell'}.
\end{equation}

The eigenvalue, $\lambda_\ell$, 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

\begin{equation}\label{eqn:sum-of-pcs}
x_{ik} = \sum_{\ell=0}^r z_{k\ell}\mkern2mu v_{i\ell}.
\end{equation}

The coefficients, $z_{k\ell}$, 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 $z_{k\ell}=\sum_{i=0}^m v_{i\ell}\mkern2mu x_{ik}$, which follows from the orthonormality of the principal components (equation <ref>).

We performed PCA separately on the frequency-averaged profile residuals in each band, using a singular value decomposition of the profile residuals. The three most significant principal components in each band are shown in Figure <ref>, along with the correspoding amplitudes, $z_{ij}$. Multiplying the amplitudes by the corresponding principal components and adding them up, as in equation (<ref>), produces a time-dependent model for the shape of the profile, which will be discussed further in Section <ref> below.

§.§ 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  [Pennucci et al., 2014, Pennucci, 2019]. First, we derived a frequency-dependent template profile 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. We then 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.5MHz 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 <ref>, and the TOA residuals are shown broken down by frequency in Figure <ref>. 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

l c c c c c c c

Degree of shape change in each band. Fractional shape differences are calculated using the root-mean-square method described in Section <ref>.
Receiver     Epoch     MJD     2cShape difference

Fixed phase     Fit phase
GBT 1500 MHz     onset     59334     $56.1\pm0.4\%$     $26.7\pm0.3\%$

GBT 820 MHz     onset     59333     $61.1\pm0.5\%$     $24.5\pm0.4\%$

CHIME 600 MHz     onset     59321–81     $42.3\pm0.6\%$     $37.1\pm0.6\%$

GBT 1500 MHz     most recent     59945     $7.08\pm0.32\%$     $2.60\pm0.32\%$

GBT 820 MHz     most recent     59916     $8.64\pm0.26\%$     $2.89\pm0.26\%$

GBT 600 MHz     most recent     59916–76     $11.3\pm1.0\%$     $10.1\pm1.0\%$

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 <ref>. 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.

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 Figure <ref>. All three portraits are de-dispersed with the same dispersion measure, 15.9638pc 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.

Overview 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 Table <ref>, 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.

Profile residuals around the shape change event, derived from the profiles in Figure <ref> by subtracting a “fixed-phase” template (cf. Section <ref>). As in Figure <ref>, 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 Figure <ref>, 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.

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

Principal component analysis of the profile residuals in each band. For each of the three bands described in Table <ref> (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.

l c c c c c c c

Time evolution of first principal component
Receiver     Model     $a$     $b$     $\tau$ (d)     $\alpha$     $\chi^2_r$
E     $1.05\pm0.06$         $271\pm12$         $512$

GBT 1500 MHz     O     $1.14\pm0.04$     $0.124\pm0.008$     $156\pm8$         $111$

P     $1.36\pm0.05$         $190\pm40$     $1.59\pm0.06$     $95.7$

E     $1.42\pm0.08$         $253\pm17$         $357$

GBT 820 MHz     O     $1.45\pm0.04$     $0.271\pm0.015$     $128\pm6$         $37.9$

P     $1.86\pm0.07$         $120\pm24$     $1.17\pm0.11$     $40.3$

E     $1.318\pm0.017$         $275\pm4$         $2.88$

CHIME 600 MHz     O     $1.308\pm0.018$     $0.206\pm0.009$     $159\pm5$         $2.00$

P     $1.580\pm0.026$         $204\pm23$     $1.46\pm0.10$     $1.95$

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 750MHz (cf. Figure <ref>), can be seen in both of these plots: the amplitude of the TOA residuals is greater in the 820MHz band than at either 600MHz or 1500MHz; and the apparent DM residuals, which roughly represent the slope of the TOA residuals as a function of frequency, are positive at 600MHz and negative at 1500MHz.

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 25MHz $\times$ 30day bins. The frequency dependence is not consistent with the $\nu^{-2}$ form expected for dispersive behavior, but instead is non-monotonic, with the TOA residuals reaching their most negative value at around 750MHz.

The average profiles seen in Figure <ref> illustrate the form of the shape change at different radio frequencies. As indicated in Table <ref>, 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. Both before and after the onset of the event, the profile shows a noticeable dependence on frequency, as can be seen in Figure <ref>. Five main emission components can be identified: in addition to the bright central peak, there are two leading components and two trailing components, one of which is evident in the original profile only as a slight change in slope on the trailing edge of the main peak, but appears much more distinct after the shape change. After the shape change, the main peak is narrower, while the leading “shoulder” component is brighter, and the trailing shoulder is broader. The second trailing component 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 <ref> demonstrates the fact that the profile has recovered toward its original shape over the months that have elapsed since the event. At 1.5GHz, 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 820MHz, 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.5GHz, 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 <ref>, which shows the profile residuals. The residuals can be seen to gradually decrease over time, while approximately maintaining their shape. As seen in Figure <ref>, this recovery can be described primarily as a change in the amplitude of the first principal component, $z_1$, as a function of 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):

\begin{equation}
z_1(t)=a e^{-(t-t_0)/\tau},
\end{equation}

a decaying exponential with persistent offset (Model O):

\begin{equation}
z_1(t) = ae^{-(t-t_0)/\tau} + b,
\end{equation}

and a power-law model (Model P):

\begin{equation}
z_1(t) = a\sbrack{1+\frac{t-t_0}{\tau}}^{-\alpha}.
\end{equation}

The results can be seen in Table <ref>. In each case, $t_0$ was taken to be MJD 59321.0 exactly; 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 <ref>. In the GBT 820MHz and 1.5GHz 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.5e-3pc cm^-3 in the 1.5GHz band and -2.5e-3pc cm^-3 in the 820MHz 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.5e-3pc 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 <ref> and <ref> 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).

§ INTERPRETATION

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 800MHz, and all three show at least some evidence for associated profile shape changes, which are only gradually frequency-dependent. Moreover, the recovery from the new event appears to have taken place at a similar rate as that of the previous events, with the profile shape and DM taking roughly a year to return to a state resembling that before the event. This is approximately twice as long as the recovery associated with each of the two previous events. It may therefore be reasonable to assume that the new event has the same physical origin as those previous events.

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–1ms), 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 $\nu^{-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 $\nu^{-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.

Many pulsars exhibit mode-changing, in which the profile abruptly changes back and forth between 2 or more shapes (e.g. Backer, 1970).  This phenomenon has been shown to extend to millisecond pulsars [Mahajan et al., 2018].  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.

The event has some features similar to those seen in 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, so the rapid onset and subsequent recovery over months seen in the shape change event in J1713+0747 are both glitch-like features. However, the time-of-arrival changes associated with glitches are independent of radio frequency, unlike those seen here, and are usually accompanied by long-term changes in pulse period or period derivative. Furthermore, glitches are not usually associated with profile shape changes. There have been some exceptions to this: 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. Additionally, observing during a glitch in the Vela pulsar, PSR B0833$-$45, Palfreyman et al., 2018 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 in each scenario is an achromatic step in the pulsar's spin frequency, which is not present in the current case. Furthermore, the shape changes previously 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, as Brook et al., 2014 argue was the case for PSR B0736$-$40 in 2006. Such an event could result in observable changes to a pulsar's radio emission if asteroidal material entered the pulsar's magnetosphere, altering the flow of particles and thus the radio emission [Cordes & Shannon, 2008]. However, this is relatively unlikely for MSPs such as PSR J1713+0747. Since the light-cylinder radius $r_{\mr{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 1500MHz, has decreased abruptly, and appears to be recovering toward its original value on a timescale of several months. Although it is frequency-dependent, 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.

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.

The NANOGrav Physics Frontiers Center receives support from the National Science Foundation (NSF) under award number 1430284. The Green Bank Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. We acknowledge that CHIME is located on the traditional, ancestral, and unceded territory of the Syilx/Okanagan people. We are grateful to the staff of the Dominion Radio Astrophysical Observatory, which is operated by the National Research Council of Canada.  CHIME is funded by a grant from the Canada Foundation for Innovation (CFI) 2012 Leading Edge Fund (Project 31170) and by contributions from the provinces of British Columbia, Québec and Ontario. The CHIME/FRB Project, which enabled development in common with the CHIME/Pulsar instrument, is funded by a grant from the CFI 2015 Innovation Fund (Project 33213) and by contributions from the provinces of British Columbia and Québec, and by the Dunlap Institute for Astronomy and Astrophysics at the University of Toronto. Additional support was provided by the Canadian Institute for Advanced Research (CIFAR), McGill University and the McGill Space Institute thanks to the Trottier Family Foundation, and the University of British Columbia. The CHIME/Pulsar instrument hardware was funded by NSERC RTI-1 grant EQPEQ 458893-2014. This research was enabled in part by support provided by WestGrid (www.westgrid.ca) and Compute Canada (www.computecanada.ca).

Support for H. T. Cromartie is provided by NASA through the NASA  Hubble  Fellowship  Program  grant  #HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. T. Dolch is supported by an NSF Astronomy and Astrophysics Grant (AAG) award number 2009468. E. C. Ferrara is supported by NASA under award number 80GSFC21M0002. M. T. Lam is supported by an NSF Astronomy and Astrophysics Grant (AAG) award number 2009468. The Flatiron Institute is supported by the Simons Foundation. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. J. W. McKee is a CITA Postdoctoral Fellow: This work was supported by Ontario Research Fund—research Excellence Program (ORF-RE) and the Natural Sciences and Engineering Research Council of Canada (NSERC) [funding reference CRD 523638-18]. T. T. Pennucci acknowledges support from the MTA-ELTE Extragalactic Astrophysics Research Group, funded by the Hungarian Academy of Sciences (Magyar Tudományos Akadémia), which was used during the development of this research. S. M. Ransom is a CIFAR Fellow. Pulsar research at UBC is supported by an NSERC Discovery Grant and by CIFAR. Portions of this work performed at NRL were supported by ONR 6.1 basic research funding.

Facilities: Green Bank Observatory, CHIME

Software: DSPSR [van Straten & Bailes, 2010, van Straten & Bailes, 2011], PSRCHIVE [Hotan et al., 2004, van Straten et al., 2011],  (<https://github.com/demorest/psrtools>),  [Demorest, 2018], Astropy [Astropy Collaboration et al., 2013, Astropy Collaboration et al., 2018, Astropy Collaboration et al., 2022], PulsePortraiture [Pennucci et al., 2016], PINT [Luo et al., 2021]

[Alam et al., 2021]
Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2021a,
, 252, 4, 10.3847/1538-4365/abc6a0

[Alam et al., 2021]
—. 2021b, , 252, 5, 10.3847/1538-4365/abc6a1

[Anderson & Itoh, 1975]
Anderson, P. W., & Itoh, N. 1975, , 256, 25, 10.1038/256025a0

[Arzoumanian et al., 2020]
Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, , 905,
L34, 10.3847/2041-8213/abd401

[Astropy Collaboration et al., 2013]
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,
, 558, A33, 10.1051/0004-6361/201322068

[Astropy Collaboration et al., 2018]
Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al.
2018, , 156, 123, 10.3847/1538-3881/aabc4f

[Astropy Collaboration et al., 2022]
Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022,
apj, 935, 167, 10.3847/1538-4357/ac7c74

[Backer, 1970]
Backer, D. C. 1970, , 228, 1297, 10.1038/2281297a0

[Backer et al., 2000]
Backer, D. C., Wong, T., & Valanju, J. 2000, , 543, 740,

[Breton et al., 2008]
Breton, R. P., Kaspi, V. M., Kramer, M., et al. 2008, Science, 321,
104, 10.1126/science.1159295

[Brook et al., 2014]
Brook, P. R., Karastergiou, A., Buchner, S., et al. 2014, , 780,
L31, 10.1088/2041-8205/780/2/L31

[Chen et al., 2021]
Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, , 508, 4970,

[CHIME Collaboration et al., 2022]
CHIME Collaboration, Amiri, M., Bandura, K., et al. 2022, , 261,
29, 10.3847/1538-4365/ac6fd9

[CHIME/Pulsar Collaboration et al., 2021]
CHIME/Pulsar Collaboration, Amiri, M., Bandura, K. M., et al. 2021,
, 255, 5, 10.3847/1538-4365/abfdcb

[Cognard & Backer, 2004]
Cognard, I., & Backer, D. C. 2004, , 612, L125,

[Cordes, 2013]
Cordes, J. M. 2013, Classical and Quantum Gravity, 30, 224002,

[Cordes & Shannon, 2008]
Cordes, J. M., & Shannon, R. M. 2008, , 682, 1152,

[Cordes et al., 2017]
Cordes, J. M., Wasserman, I., Hessels, J. W. T., et al. 2017, ,
842, 35, 10.3847/1538-4357/aa74da

[Cornish & Sesana, 2013]
Cornish, N. J., & Sesana, A. 2013, Classical and Quantum Gravity, 30,
224005, 10.1088/0264-9381/30/22/224005

[Demorest, 2007]
Demorest, P. B. 2007, PhD thesis, University of California, Berkeley.
<https://www.cv.nrao.edu/ pdemores/thesis.pdf>

[Demorest, 2018]
—. 2018, nanopipe: Calibration and data reduction pipeline for pulsar
timing, Astrophysics Source Code Library, record ascl:1803.004.

[Demorest et al., 2013]
Demorest, P. B., Ferdman, R. D., Gonzalez, M. E., et al. 2013, ,
762, 94, 10.1088/0004-637X/762/2/94

[Desvignes et al., 2016]
Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, , 458,
3341, 10.1093/mnras/stw483

[Dolch et al., 2014]
Dolch, T., Lam, M. T., Cordes, J., et al. 2014, , 794, 21,

[DuPlain et al., 2008]
DuPlain, R., Ransom, S., Demorest, P., et al. 2008, in Society of
Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7019,
Advanced Software and Control for Astronomy II, ed. A. Bridger & N. M.
Radziwill, 70191D, 10.1117/12.790003

[Goncharov et al., 2021]
Goncharov, B., Reardon, D. J., Shannon, R. M., et al. 2021, ,
502, 478, 10.1093/mnras/staa3411

[Graham Smith et al., 2011]
Graham Smith, F., Lyne, A. G., & Jordan, C. 2011, , 410, 499,

[Haskell & Melatos, 2015]
Haskell, B., & Melatos, A. 2015, International Journal of Modern Physics
D, 24, 1530008, 10.1142/S0218271815300086

[Hellings & Downs, 1983]
Hellings, R. W., & Downs, G. S. 1983, , 265, L39,

[Hobbs et al., 2019]
Hobbs, G., Dai, S., Manchester, R. N., et al. 2019, Research in
Astronomy and Astrophysics, 19, 020, 10.1088/1674-4527/19/2/20

[Hotan et al., 2004]
Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, , 21, 302,

[Jennings et al., 2022]
Jennings, R. J., NANOGrav Collaboration, & CHIME/Pulsar Collaboration.
2022, The Astronomer's Telegram, 15223, 1.

[Joshi et al., 2018]
Joshi, B. C., Arumugasamy, P., Bagchi, M., et al. 2018, Journal of
Astrophysics and Astronomy, 39, 51, 10.1007/s12036-018-9549-y

[Karastergiou et al., 2011]
Karastergiou, A., Roberts, S. J., Johnston, S., et al. 2011, ,
415, 251, 10.1111/j.1365-2966.2011.18697.x

[Keith et al., 2013]
Keith, M. J., Shannon, R. M., & Johnston, S. 2013, , 432, 3080,

[Kerr, 2015]
Kerr, M. 2015, , 452, 607, 10.1093/mnras/stv1296

[Kerr et al., 2020]
Kerr, M., Reardon, D. J., Hobbs, G., et al. 2020, , 37, e020,

[Kramer, 1998]
Kramer, M. 1998, , 509, 856, 10.1086/306535

[Kramer et al., 2006]
Kramer, M., Lyne, A. G., O'Brien, J. T., Jordan, C. A., & Lorimer, D. R. 2006,
Science, 312, 549, 10.1126/science.1124060

[Lam, 2021]
Lam, M. T. 2021, Research Notes of the American Astronomical Society, 5, 167,

[Lam et al., 2016]
Lam, M. T., Cordes, J. M., Chatterjee, S., et al. 2016, , 821, 66,

[Lam et al., 2018]
Lam, M. T., Ellis, J. A., Grillo, G., et al. 2018, , 861, 132,

[Lentati et al., 2015]
Lentati, L., Alexander, P., & Hobson, M. P. 2015, , 447, 2159,

[Lentati et al., 2017]
Lentati, L., Kerr, M., Dai, S., et al. 2017, , 466, 3706,

[Levin et al., 2016]
Levin, L., McLaughlin, M. A., Jones, G., et al. 2016, , 818, 166,

[Lin et al., 2021]
Lin, F. X., Lin, H.-H., Luo, J., et al. 2021, , 508, 1115,

[Luo et al., 2021]
Luo, J., Ransom, S., Demorest, P., et al. 2021, , 911, 45,

[Lyne et al., 2010]
Lyne, A., Hobbs, G., Kramer, M., Stairs, I., & Stappers, B. 2010, Science,
329, 408, 10.1126/science.1186683

[Mahajan et al., 2018]
Mahajan, N., van Kerkwijk, M. H., Main, R., & Pen, U.-L. 2018, ,
867, L2, 10.3847/2041-8213/aae713

[McKee et al., 2016]
McKee, J. W., Janssen, G. H., Stappers, B. W., et al. 2016, ,
461, 2809, 10.1093/mnras/stw1442

[Meyers & CHIME/Pulsar Collaboration, 2021]
Meyers, B., & CHIME/Pulsar Collaboration. 2021, The Astronomer's Telegram,
14652, 1.

[Nan & Li, 2013]
Nan, R., & Li, D. 2013, in Materials Science and Engineering Conference
Series, Vol. 44, Materials Science and Engineering Conference Series, 012022,

[Nobleson et al., 2022]
Nobleson, K., Agarwal, N., Girgaonkar, R., et al. 2022, , 512,
1234, 10.1093/mnras/stac532

[Noutsos et al., 2020]
Noutsos, A., Desvignes, G., Kramer, M., et al. 2020, , 643, A143,

[Osłowski et al., 2011]
Osłowski, S., van Straten, W., Hobbs, G. B., Bailes, M., &
Demorest, P. 2011, , 418, 1258,

[Palfreyman et al., 2018]
Palfreyman, J., Dickey, J. M., Hotan, A., Ellingsen, S., & van
Straten, W. 2018, , 556, 219, 10.1038/s41586-018-0001-x

[Pennucci, 2019]
Pennucci, T. T. 2019, , 871, 34, 10.3847/1538-4357/aaf6ef

[Pennucci et al., 2014]
Pennucci, T. T., Demorest, P. B., & Ransom, S. M. 2014, , 790, 93,

[Pennucci et al., 2016]
—. 2016, Pulse Portraiture: Pulsar timing, Astrophysics Source Code
Library, record ascl:1606.013.

[Perera et al., 2019]
Perera, B. B. P., DeCesar, M. E., Demorest, P. B., et al. 2019, ,
490, 4666, 10.1093/mnras/stz2857

[Prestage et al., 2015]
Prestage, R. M., Bloss, M., Brandt, J., et al. 2015, in 2015 URSI-USNC
Radio Science Meeting, 4, 10.1109/USNC-URSI.2015.7303578

[Reardon et al., 2021]
Reardon, D. J., Shannon, R. M., Cameron, A. D., et al. 2021, ,
507, 2137, 10.1093/mnras/stab1990

[Roshi et al., 2011]
Roshi, D. A., Bloss, M., Brandt, P., et al. 2011, in 2011 XXXth URSI General
Assembly and Scientific Symposium, 1–4,

[Shannon et al., 2016]
Shannon, R. M., Lentati, L. T., Kerr, M., et al. 2016, , 828, L1,

[Shaw et al., 2022]
Shaw, B., Stappers, B. W., Weltevrede, P., et al. 2022, , 513,
5861, 10.1093/mnras/stac1156

[Singha et al., 2021]
Singha, J., Surnis, M. P., Joshi, B. C., et al. 2021, , 507, L57,

[Spiewak et al., 2022]
Spiewak, R., Bailes, M., Miles, M. T., et al. 2022a,
, 39, e027, 10.1017/pasa.2022.19

[Spiewak et al., 2022]
—. 2022b, arXiv e-prints, arXiv:2204.04115.

[Splaver et al., 2005]
Splaver, E. M., Nice, D. J., Stairs, I. H., Lommen, A. N., & Backer,
D. C. 2005, , 620, 405, 10.1086/426804

[Stairs et al., 2019]
Stairs, I. H., Lyne, A. G., Kramer, M., et al. 2019, , 485, 3230,

[Taylor, 1992]
Taylor, J. H. 1992, Philosophical Transactions of the Royal Society of London
Series A, 341, 117, 10.1098/rsta.1992.0088

[Turner et al., 2021]
Turner, J. E., McLaughlin, M. A., Cordes, J. M., et al. 2021, ,
917, 10, 10.3847/1538-4357/abfafe

[van Straten & Bailes, 2010]
van Straten, W., & Bailes, M. 2010, DSPSR: Digital Signal Processing
Software for Pulsar Astronomy, Astrophysics Source Code Library, record

[van Straten & Bailes, 2011]
—. 2011, , 28, 1, 10.1071/AS10021

[van Straten et al., 2011]
van Straten, W., Demorest, P., Khoo, J., et al. 2011, PSRCHIVE:
Development Library for the Analysis of Pulsar Astronomical Data,
Astrophysics Source Code Library, record ascl:1105.014.

[Weisberg et al., 1989]
Weisberg, J. M., Romani, R. W., & Taylor, J. H. 1989, , 347, 1030,

[Weltevrede et al., 2011]
Weltevrede, P., Johnston, S., & Espinoza, C. M. 2011, , 411, 1917,

[Xu et al., 2021]
Xu, H., Huang, Y. X., Burgay, M., et al. 2021, The Astronomer's
Telegram, 14642, 1.

[Zhu et al., 2019]
Zhu, W. W., Desvignes, G., Wex, N., et al. 2019, , 482, 3249,