1. Introduction
Millisecond pulsars (MSPs – Backer et al. Reference Backer, Kulkarni, Heiles, Davis and Goss1982; Alpar et al. Reference Alpar, Cheng, Ruderman and Shaham1982) are extraordinarily stable rotators (Verbiest et al. Reference Verbiest2009), comparable to terrestrial atomic clocks over years to decades (Taylor Reference Taylor1991). While a large fraction of canonical pulsars exhibit large rotational irregularities (Hobbs, Lyne, & Kramer Reference Hobbs, Lyne and Kramer2010; Lower et al. Reference Lower2025) or frequent glitching (Basu et al. Reference Basu2022), these phenomena are much more rare, or lower in magnitude, in MSPs (Shannon & Cordes Reference Shannon and Cordes2010; McKee et al. Reference McKee2016).
Time-integrated MSP pulse profiles have in general proven to be exceptionally stable compared with the canonical pulsar population (Liu et al. Reference Liu2012), reflecting the stability of their magnetospheres. The rotational and profile stability of MSPs can be exploited to carry out high-precision pulsar timing experiments, including precise pulsar mass measurements (Choudhury et al. Reference Choudhury2024; Reardon et al. Reference Reardon2024; Rutherford et al. Reference Rutherford2024), stringent tests of general relativity (Kramer et al. Reference Kramer2021a,Reference Kramerb) and searches for gravitational waves (Agazie et al. Reference Agazie2023a; EPTA Collaboration et al. Reference Collaboration2023b; Reardon et al. Reference Reardon2023; Xu et al. Reference Xu2023; Miles et al. Reference Miles2025b).
Despite the apparent stability of MSP integrated profiles, some MSPs have exhibited similar temporal profile variations to those seen in canonical pulsars. One class of pulsar profile variability is ‘jitter’ (Parthasarathy et al. 2021), which refers to variations of pulsar profile shapes on the timescale of individual rotations. This may lead to small variations in integrated pulse profiles, particularly for short observations. Short-term discrete profile changes, such as mode-changing and nulling (Miles et al. Reference Miles2022; Nathan et al. Reference Nathan2023), can also occur on short timescales, ranging from individual rotations to days. Stochastic profile variability on longer timescales of months – years has been observed in both canonical pulsars (Brook et al. Reference Brook2016; Shaw et al. Reference Shaw2022; Lower et al. Reference Lower2023; Basu et al. Reference Basu2024; Lower et al. Reference Lower2025) and MSPs (Backer & Sallmen Reference Backer and Sallmen1997; Osłowski et al. Reference Osłowski, van Straten, Hobbs, Bailes and Demorest2011; Brook et al. Reference Brook2018; Padmanabh et al. Reference Padmanabh2021; Nathan et al. Reference Nathan2023; Fiore et al. Reference Fiore2025), and is thought to arise due to long-term variability in the pulsar magnetosphere. Additionally, propagation effects through the interstellar medium, such as scattering (Lewandowski, Kowalińska, & Kijak Reference Lewandowski, Kowalińska and Kijak2015), diffractive scintillation (Liu et al. Reference Liu2022), and plasma lensing (Graham Smith, Lyne, & Jordan Reference Graham Smith, Lyne and Jordan2011), can cause profile variability, although these are clearly not intrinsic to the pulsar itself.
Discrete profile change events – the subject of this paper – manifest as a sudden augmentation of pulse profile components or the appearance of new components that decay over time. There are several noteworthy examples of long-term discrete profile changes reported in the MSP population: PSR J0437
$-$
4715 (Goncharov et al. Reference Goncharov2021), PSR J1643
$-$
1224 (Shannon et al. Reference Shannon2016; Goncharov et al. Reference Goncharov2021), and two prior events on PSR J1713
$+$
0747 (hereafter ‘J1713’) which we describe below. Transient profile changes have also been observed from several canonical pulsars (Lyne et al. Reference Lyne, Hobbs, Kramer, Stairs and Stappers2010; Karastergiou et al. Reference Karastergiou2011; Brook et al. Reference Brook2014; Brook et al. Reference Brook2016; Shaw et al. Reference Shaw2022; Keith et al. Reference Keith2025). However, most of these transient profile events are different in character to the abrupt change followed by gradual decays seen in MSPs like J1713 – for example, they may arise more gradually, manifest as one component of a long-term ‘state-switching’ phenomenon, or do not exhibit a gradual decay, instead forming part of a continuous process of profile variability. Close analogues to the abrupt profile change event like PSR J1713
$+$
0747 have been observed in canonical pulsars like PSRs J0738
$-$
4042 (Karastergiou et al. Reference Karastergiou2011; Brook et al. Reference Brook2014), B2035
$+$
36 (Shaw et al. Reference Shaw2022; Keith et al. Reference Keith2025), and J1119
$-$
6127, which showed erratic pulse shape changes shortly after large glitches (Weltevrede, Johnston, & Espinoza Reference Weltevrede, Johnston and Espinoza2011) and X-ray outbursts (Dai et al. Reference Dai2018). Given the apparent long-term stability of most MSP magnetospheres, studies of this class of profile change event in MSPs may yield important clues into how the phenomenon occurs across the population as a whole, and to the emission mechanism itself.
PSR J1713
$+$
0747 (Foster, Wolszczan, & Camilo Reference Foster, Wolszczan and Camilo1993) is a 4.57-millisecond-period MSP in a binary system with a white dwarf companion. Its long-term timing stability, bright, narrow pulse profile, and equatorial sky position has made it a centrepiece of many PTA observing programmes (Agazie et al. Reference Agazie2023b; EPTA Collaboration et al. Reference Collaboration2023a; Zic et al. Reference Zic2023; Xu et al. Reference Xu2023; Miles et al. Reference Miles2025a).
J1713 has undergone two frequency-dependent fluctuations in its times of arrival which have been attributed to profile change events. The first, in October 2008, lasted between 100 and 200 days before recovering (Coles et al. Reference Coles2015; Zhu et al. Reference Zhu2015; Desvignes et al. Reference Desvignes2016). The second event occurred
$\sim$
8 yr later, in April 2016 (Lam et al. Reference Lam2018). Lin et al. (Reference Lin2021) suggested both events likely originated from discrete dispersion measure (DM) variations. On the other hand, Goncharov et al. (Reference Goncharov2021) concluded that the 2016 event was likely magnetospheric in origin. Although there is no reported evidence for pulse shape changes during the 2008 event, Jennings et al. (Reference Jennings2024) proposed that both historical events were likely related to the magnetosphere.
In mid-April 2021, J1713 exhibited a significant change to its radio pulse profile (Xu et al. Reference Xu2021; Meyers & Chime/Pulsar Collaboration 2021; Singha et al. Reference Singha2021b). The time of the profile change was determined to be between Modified Julian Day (MJD) 59320 – 59321 (Meyers & Chime/Pulsar Collaboration 2021). Studies of the April 2021 event by Singha et al. (Reference Singha2021a), Jennings et al. (Reference Jennings2024) and Wang et al. (Reference Wang2024) have revealed complex, frequency-dependent profile shape changes, suggesting a magnetospheric origin as opposed to propagation effects. These studies observed the profile recovery over time, reporting different timelines for its return to a near pre-event state (Singha et al. Reference Singha2021a; Wang et al. Reference Wang2024; Jennings et al. Reference Jennings2024). Both Singha et al. (Reference Singha2021a) and Wang et al. (Reference Wang2024) also reported that higher frequencies exhibited more pronounced changes. These studies collectively disfavour alternative explanations such as glitches, DM changes, scintillation and jitter (Lam Reference Lam2021; Jennings et al. Reference Jennings2024).
Radio emission from pulsars is highly polarised, exhibiting frequency-dependent linear and circular polarisation profiles that have distinct morphological features (Lyne & Smith Reference Lyne and Smith1968; Gould & Lyne Reference Gould and Lyne1998). The polarisation profiles encode crucial information about the emission mechanism and propagation of radiation through the pulsar magnetosphere (Kramer et al. Reference Kramer1999; Yan et al. Reference Yan2011; Karastergiou et al. Reference Karastergiou2024; Xu et al. Reference Xu2025). Linear polarisation is particularly important, as the observed polarisation position angle (PA) can be used to infer the magnetic field geometry of the pulsar via the ‘rotating vector model’ (Radhakrishnan et al. Reference Radhakrishnan, Cooke, Komesaroff and Morris1969). While many pulsars show smooth PA tracks, some pulsars can also exhibit abrupt 90* jumps across pulse phase, corresponding to transitions between the orthogonal polarisation modes (OPMs; Manchester Reference Manchester1975).
To date, no detailed study of the polarimetric behaviour of the J1713 profile change has been conducted. Discrete changes in polarisation properties, such as those induced by profile change events, therefore offer important insights into evolving magnetospheric conditions and sources of instability that may affect high-precision timing experiments. Wideband spectro-polarimetric observations are particularly crucial for probing these phenomena, offering access to information about the magnetosphere that would otherwise be inaccessible (Oswald, Karastergiou, & Johnston Reference Oswald, Karastergiou and Johnston2020; Oswald et al. Reference Oswald2023b; Lower et al. Reference Lower2024). Furthermore, because J1713 is one of the most important MSPs in global PTA campaigns, detailed characterisation of this profile change event is critical for informing mitigation strategies in precision timing analysis such as searches for low-frequency gravitational waves.
In this paper, we present a wideband polarimetric analysis of the ongoing recovery of the profile of J1713 following the April 2021 change event – the most significant profile change yet recorded on an MSP. Under the Parkes Pulsar Timing Array observing programme, J1713 has been regularly observed using the Ultra-Wideband Low-frequency receiver (UWL; Hobbs et al. Reference Hobbs2020) on Murriyang, CSIRO’s Parkes Radio Telescope. With a frequency band spanning 704–4 032 MHz, well-calibrated polarimetry, and a long-term continuous wide-band dataset spanning
$\gt 3$
yr, this dataset offers the opportunity to characterise this profile change in greater depth than previous studies. Observation details and methods are presented in Section 2, analysis and findings in Section 3, and further discussions in Section 4.
2. Observations and methods
2.1 Observations and data reduction
Our observations were taken under two separate observing programmes on Murriyang: The Parkes Pulsar Timing Array programme (PPTA; Manchester et al. Reference Manchester2013, project code P456), and a separate, targeted programme (P1172). The PPTA observes a set of 37 MSPs with approximately 3-week cadence (see Zic et al. Reference Zic2023 for details of the observing programme). Under the P1172 programme, we took targeted high-cadence (
$\sim$
2–7 days) observations of J1713 to more closely monitor the recovery from the April 2021 profile change event, between October 2022 and April 2024. All observations were carried out using the UWL receiver, recorded in pulsar fold mode and coherently de-dispersed to the pulsar DM value
$15.988$
pc cm
$^{-3}$
. The observations were typically 64 min in duration, although occasionally observations were cut short for operational or instrumental reasons.
The observations were flagged, calibrated, and reduced following the same procedures used for the PPTA Third Data Release (PPTA DR3), as described in Zic et al. (Reference Zic2023). In this work, we make use of the sub-banded data products, using the same frequency sub-band divisions defined in the PPTA DR3, which we list in Table 1 for clarity. We installed the timing ephemeris from the PPTA DR3 onto all observation archive files. To ensure sufficient signal-to-noise ratio (S/N), we frequency-averaged the profiles within each of the 8 sub-bands. For each sub-band, we removed any observations with a S/N less than 50. We also excluded data with corrupted polarisation measurements due to digitiser synchronisation errors; this affected small frequency ranges for short amounts of time, most significantly around December 2021–February 2022 across 1.0–1.3 GHz and February 2022–May 2022 from 3.0–3.5 GHz.
UWL sub-band ranges (where
$\nu_{\mathrm{min}}$
and
$\nu_{\mathrm{max}}$
represent the minimum and maximum frequency), bandwidth and centre frequency (
$\nu_{c})$
.
In this work, we analyse changes to linear polarisation PAs over time. This requires measurement and correction for variations in the interstellar and ionospheric rotation measure (RM), in addition to the data processing performed in the PPTA DR3. To this end, we used the rmfit tool in the psrchive software suite (Hotan, van Straten, & Manchester Reference Hotan, van Straten and Manchester2004; van Straten & Bailes Reference van Straten and Bailes2011) to measure the RM for each observation, manually setting the RM bounds from -10 to 30
$\text{rad m}^{-2}$
, and using 40 steps. We then used the general-purpose pulsar data manipulation tool pam to correct the profile archives for RM at an infinite reference frequency, using the --aux_rm option.
2.2 Profile alignment
The dramatic and sustained change in the profile shape of J1713 (see Figures 1, 2) means that one of the fundamental assumptions of pulsar timing – that the pulse shape remains stable over time – is violated. This leads to a degeneracy between profile variability (e.g. due to magnetospheric changes) and pulse arrival times, because a mismatch between a reference template and the true profile can be erroneously modeled as a profile phase shift. Conversely, if there were variations to the rotational properties of the pulsar, the change in profile shape makes it challenging to disentangle rotational effects from those that purely augment the observed pulse profile (e.g. in the pulsar magnetosphere). For the pulse profile analyses presented in this paper, one possible assumption is that the intrinsic timing properties of the pulsar have not changed. In this case, the expected rotational phase location of the radio pulse is determined by the ephemeris prediction. While this is subject to usual variations that affect pulsar timing (e.g. DM variations, timing noise, secular changes in orbital parameters, etc.), using an established timing ephemeris to predict rotational phase allows careful comparison of specific sub-components of the pulse profile which may (or may not) have changed significantly. On the other hand, if we relax the assumption that the intrinsic rotational properties did not change, the profile change can be studied by minimising the difference between the observed pulse profile and a standard template. This necessitates aligning the observed profiles with this standard template, and, by design, minimises the difference between the observed profile and the template. In practice, we achieve this by using the Taylor (Reference Taylor1992) Fourier phase-gradient alignment algorithm, fitting the phase shift of the profiles using the templates which we describe in Section 2.3. A genuine change to the intrinsic rotation wouldn’t produce the profile residuals that we observe in the FFT-aligned methodoly we have applied, and the appearance of such residuals are the result of intrinsic profile shape changes. In the analyses presented here, we employ both the ephemeris-predicted, and template-aligned methods to predict the rotational phase location of the pulse profiles over time.
Stokes I, PA, Linear polarisation and Stokes V profiles colour-mapped as a function of frequency per row. The first column represents the template profiles, the second column represents MJD 59368 (47 days post-event), and the last column represents MJD 60049 (roughly two years post-event). Annotated are five profile components: (A) the leading peak, (B) the profile shoulder, (C) the main pulse peak, (D) the descending gradient, and (E) the trailing peak. The four OPM sub-components are also annotated. In addition to the normalisation by the integrated Stokes I flux density, we additionally normalise the Stokes I, L, and Stokes V profiles for each epoch in this plot by the maximum of normalised Stokes I flux densities across all sub-bands, scaling the normalised intensities to a maximum of 1.
Stokes I, L and Stokes V profile residuals as a function of time, where pulse phase at each epoch is predicted from the timing ephemeris. The top row shows Stokes I, the middle row shows linearly polarised intensity (“L”), and the bottom row shows Stokes V. The red dashed line indicates the time of the profile change event. The colour bar indicates the intensity of the profile residual. As can be seen, Stokes I and L were affected, whereas Stokes V was minimally affected in the lower sub-bands A–E, but shows variations in the higher-frequency sub-bands sbF – sbH. Additionally, the profile has not yet returned to its pre-event state, indicating either a long recovery timescale or potential reconfiguration of the profile.
2.3 Sub-banded template generation
We used all available UWL observations prior to the event date to construct frequency-averaged templates of all four polarisations (I, Q, U, V) for each of the constructed sub-bands. These observations spanned from October 2018 to the end of March 2021.
The polarised, sub-band profile templates were generated as follows. We first estimate the total intensity S/N in every pre-event observation using psrstat and sort these profiles in descending order according to their S/N. We then begin an iterative process of constructing the templates, beginning with the highest S/N profile as the initial reference. We use the Taylor (Reference Taylor1992) Fourier phase-gradient alignment algorithm to align the second-highest S/N profile to the first. This phase shift is computed on the total intensity profile but applied to all four Stokes parameters (I, Q, U, V). We then construct a new template by taking the weighted mean of all of the aligned profiles, which is then used as the template profile in the next step. We apply this process iteratively until all profiles are aligned. The template thus generated is the mean of all aligned profiles.
This template is used to produce profile residuals both for the ephemeris-predicted and template-aligned profiles. For the ephemeris-predicted dataset, we adjust the absolute phase of the template by taking a simple un-aligned mean of all pre-event profiles, calculating the phase offset between our constructed template and this simple mean, and shifting the template by this offset. This ensures that the profile residuals on the ephemeris-predicted profiles are not biased by any timing deviations present in the highest S/N profile.
We construct the template and sample profiles for the total linearly polarised intensity (hereafter L) using the polarised intensity image-based debiasing technique detailed in Müller et al. (Reference Müller, Beck and Krause2017), modified to operate on our one-dimensional profiles. This is used to mitigate the Ricean noise statistics that bias naïve calculations of
$L = \sqrt{Q^2 + U^2}$
when using noisy data. This method provides superior performance at low profile S/N compared to the method of Everett & Weisberg (Reference Everett and Weisberg2001) commonly used in the pulsar literature, which truncates L profiles below a S/N of 1.57. The Müller et al. (Reference Müller, Beck and Krause2017) debiasing method provides a continuous estimate of L regardless of the profile S/N, while also preserving the zero-mean gaussian noise statistics of the input Stokes Q and U profiles in the final L profile. The linear polarisation PA of the template, and sample profiles are computed as
$\text{PA} = 1/2~\text{arctan}2(U,Q)$
, where
$\text{arctan}2$
is the two-argument arctangent.
We note that accurately modelling the frequency evolution of MSP profiles is challenging, and we do not specifically attempt to do so in this study. However, the stability of the pre-event residuals across all sub-bands over 2.4 yr prior to the event (see Figure 2), demonstrates that our sub-banded templates, constructed from long-term average profiles, are reliable and fit-for-purpose.
2.4 Normalisation and construction of profile residuals
Processes such as scintillation and intrinsic flux density modulation can induce frequency-dependent variability in pulsar profiles, including that of J1713 (Bogdanov et al. Reference Bogdanov, Pruszyńska, Lewandowski and Wolszczan2002; Brook et al. Reference Brook2018). To compare pulse profiles across time, we therefore first normalise the template and observed profiles using the integrated flux density under the on-pulse region of its Stokes I profile. We use a fixed on-pulse window over pulse phase bins 0.38–0.68 to mitigate potential bias arising from the S/N dependence in determining the on-pulse window. The normalisation procedure ensures that variations in the profile shape, rather than its amplitude, become the primary focus of our analysis.
In some epochs, we found a small flux density offset between the peak of the normalised template and normalised observation profile, particularly in low S/N epochs, which introduced a small artificial residual in the pre-event data. We corrected this by performing a least-squares fit to determine the optimal scale factor between the aligned profile and template within the central pulse component (phase 0.462–0.562). These scale factors were then stored and also applied to the ephemeris-predicted profiles, ensuring consistent flux density scaling and minimising normalisation-induced artifacts in the residuals.
We constructed the profile residuals by subtracting the Stokes I, Stokes V, or L template from the corresponding observation. We carefully validated the data by inspecting each profile and its residuals by eye, and produced histograms of the off-pulse regions of the profile to identify observations corrupted by excessive interference or instrumental issues. Any suspect data were inspected, then accepted, reprocessed or rejected as applicable.
3. Results
In this section, we present our results, starting with an analysis of profile shape changes for the different Stokes parameters, followed by a principal component analysis (PCA) of the profile residuals. We then examine variability in the linear polarisation, PA, and orthogonal polarisation modes.
3.1 Profile shape variability
In Figure 1, we present the template for J1713 (first column) as well as two sample high S/N profiles that highlight the profile change effects at two epochs – MJD 59368 which is 47 days after the event, and MJD 60049, approximately two years after the event. Also included are the PA data. The frequency across the UWL spanning the eight sub-bands is colour-mapped with lower frequencies as the darker curves and higher frequencies as lighter curves, and Stokes I, PA, L, and Stokes V are plotted on separate rows. Annotated are five profile components, labelled: (A) the leading peak, (B) the profile shoulder, (C) the main pulse peak, (D) the descending gradient, and (E) the trailing peak. At early times after the profile change event, the central peak of the Stokes I profile (profile component C) significantly narrows across the entire band, whilst profile component D broadens into a plateau across pulse phase – increasing the width of the profile base and forming a distinct horizontal curve before resuming its decline. This horizontal feature is visible across all sub-bands, except sbG and sbH. Above approximately 1 500 MHz, the leading shoulder (profile component B at pulse phase bins
$\sim$
0.45–0.50) intensifies with increasing frequency.
The linear polarisation template exhibits two main peaks close to the centre, with two smaller trailing peaks. Each of these structures corresponds to a different OPM sub-component (see Section 3.2 for further details). These different OPM sub-components have been labelled in the first PA sub-plot. As Figure 1 highlights, there were dramatic, complex frequency-dependent changes in the relative heights of the linear polarisation profile components. These correlate with changes in the PA – with sub-component 2.1 narrowing, and sub-component 1.2 forming modal bridges to its neighbouring component. Further details of these variations are discussed in Section 3.2.2. A prominent change in the first peak of the linear polarisation profile is also observed, where the peak height increases dramatically in a non-monotonic frequency-dependent manner. As shown in the central sub-plot for linear polarisation of Figure Section 1, which represents the epoch soon after the event, the higher-frequency bands exhibit an increase by a factor of
$\sim$
2.5 compared to pre-event levels, before gradually returning to baseline values. This first peak corresponds to the OPM sub-component 1.1 (see Section 3.2.1 for further details).
While minor variations are seen in Stokes V in this figure, the overall flux density and scale of change remain much smaller than Stokes I and linear polarisation.
Figures 2 and 3 show the evolution of the profile residuals over time, spanning from the beginning of our UWL observations in October 2018 (MJD 58426) and extending to the start of June 2024 (MJD 60462). Figure 2 shows the profile residuals from ephemeris-predicted profiles, and Figure 3 shows the template-aligned profiles. In both sets of residuals, we observe that Stokes I and linear polarisation exhibit significant changes compared to the pre-event data across the full band. In contrast, Stokes V, which has the lowest flux density and therefore lowest S/N, does not show a significant deviation post-event at frequencies below
$\sim$
1 400 MHz (sub-bands A – E), with fewer than 4% of post-event observations having peak residual values above
$3\sigma$
significance. On the other hand, inspection of the post-event residuals reveals significant profile variations at higher frequencies, as can be seen in Figures 2 and 3, with 13–21% of observations showing residual peaks above
$3\sigma$
significance, across sbF – sbH (1 376–3 312 MHz).
As in Figure 2, but for the case where the profiles have been aligned to the template.
Additionally, new profile features emerge across the entire band. While most of these features gradually diminish over time, some – particularly in the higher frequency bands (
$\nu_{c}$
: 2 320 MHz; 3 312 MHz) – persist. Notably, a stable positive residual feature at approximately pulse phase 0.42 in the highest sub-band remains present through to our final observation, indicating the emergence of a new structure in the pulse profile that persists over long timescales. The overall changes in the linear polarisation profile residuals are narrower in pulse phase compared to those in Stokes I, yet still span the entire UWL bandwidth.
The residuals from our template-aligned profiles (Figure 3) show consistent results with the residuals from ephemeris-predicted profiles. New profile features are also observed under this alignment method across the band, gradually diminishing over time. The positive feature at pulse phase 0.42 in the highest sub-bands remains evident in this dataset. Using this alignment method, we also observe a sudden right-ward jump in phase followed by a gradual left-ward drift in the central components for approximately one year following the event, which is not observed in the ephemeris-aligned residuals. This is a non-physical artifact resulting from the template alignment process, as this method is sensitive to time-dependent shape changes. Due to this introduced bias, we consider the template-alignment method sub-optimal relative to the ephemeris-aligned method.
3.1.1 Principal component analysis on profile variability
Following Jennings et al. (Reference Jennings2024), we performed a PCA on both the Stokes I and linear polarisation residuals (using the sklearn library, Pedregosa et al. Reference Pedregosa2011), in order to quantify and model the ongoing recovery of the profile. In essence, this technique quantifies the complex pulse phase-dependent variations using a small number of PCs. It works by decomposing the covariance matrix of the profile residuals into orthogonal PCs (eigenvectors) that are ranked by the degree of the residual structure they explain (i.e. ranked by the magnitude of the corresponding eigenvalues). We conducted PCA on each individual sub-band, yielding PC scores that quantify the magnitude of each of the PCs over time.
We assessed the first ten PCs for each sub-band, finding that only the top three were significant. Of these, the first PC dominated, with an explained variance of 58% to 85% in Stokes I and 53% to 77% for linear polarisation among the top three PCs. Sub-bands A–C exhibited lower explained variances in both Stokes I and linear polarisation, while sub-bands D–H showed higher explained variances, which we hypothesise is due to the influence of residual radio frequency interference (RFI) remaining in the low-frequency data. We show the first PC scores (PC1 scores) across frequency and time in Figure 4. We find that following the event, the PC1 scores show an abrupt increase before gradually decaying. Motivated by this behaviour, we assessed an exponential and a power-law decay model that phenomenologically describe the long-term behaviour of PC1 scores.
Results of PCA analysis showing the first principal component scores in Stokes I(red) and linear polarisation (blue) as a function of time. The top subpanel of each plot shows the fit to the PC score data, and the bottom subpanel shows the fit residuals. A solid line indicates that the power-law model is preferred, where a dashed line indicates the exponential model. The orange dashed line indicates the profile change event, and the surrounding shaded grey region indicates the period of no observations. The
$\chi^{2}_{\mathrm{r}}$
for both models are also shown in each sub-plot.
The exponential decay model is given by
\begin{equation}D_e(t) =\begin{cases} B_{\mathrm{pre}} & \text{if } t \lt T_e \\ A_e \exp\left(-(t - T_{e})/\tau_e\right)+ B_{\mathrm{post}}, & \text{if } t \gt T_e \end{cases},\end{equation}
where
$B_{\mathrm{pre}}$
(
$B_{\mathrm{post}}$
) is the pre- (post-)event baseline level, A is the exponential amplitude,
$T_{e}$
is the event time (MJD 59320; fixed),
$\tau$
is the characteristic timescale of the exponential decay. We chose to model the pre- and post-event baseline levels separately to assess the evidence for a permanent reconfiguration of the profile, which would correspond with
$B_{\mathrm{pre}} \neq B_{\mathrm{post}}$
.
The power-law decay model is
\begin{equation}D_{pl}(t) =\begin{cases}B_{\mathrm{pre}}, & \text{for } t \lt T_e \\A_{pl} \left[1 + \dfrac{t - T_e}{\tau_{pl}} \right]^{-\alpha} + B_{\mathrm{post}}, & \text{for } t \geq T_e\end{cases}\end{equation}
where
$A_{pl}$
is the amplitude of the power-law decay at the event time,
$\alpha$
is the power-law index,
$T_e = 59\,320$
is the MJD of the profile change event,
$\tau_{pl} = 365.25$
is the reference timescale (set to one year), and
$B_{\mathrm{pre}}$
is the baseline offset.
We performed non-linear least-squares optimisation on the exponential decay and power-law models on the PC1 scores using Scipy’s curve_fit routine, with the Trust Region Reflective algorithm. We set parameter bounds to constrain the fits for both models. In the power law model, we set the bounds of A and
$\alpha$
to be positive to avoid non-physical solutions. We also did the same in the exponential model for the A and
$\tau$
parameters. We allowed both negative and positive values for both
$B_{pre}$
and
$B_{post}$
. We performed the fits for both Stokes I and linear polarisation. We estimated the uncertainties on the data points using the standard deviation of the pre-event PC1 scores.
To assess whether the exponential (exp.) or power-law (PL) model better describes the long-term evolution of the profile change, we calculated the reduced chi-squared (
$\chi^{2}_{\mathrm{r}}$
) statistic to assess overall goodness of fit and the Bayesian Information Criterion (BIC) to enable a direct model comparison between the exp. and PL models. These metrics were calculated for each sub-band and applied to both Stokes I and linear polarisation.
The BIC values are calculated as:
\begin{equation}\mathrm{BIC} = k \ln(n) + n \ln\left(\frac{\sum_{i=1}^{n} (P_{i} - D(t_{i}))^2}{n}\right)\end{equation}
where k is the number of free parameters in the model, n is the number of epochs, and
$P_{i}$
is the PC score at epoch
$t_i$
,
$D(t_{i})$
is the corresponding model prediction, and the difference of these is the residuals.
We evaluated these metrics using two versions of the PCA timeseries dataset: the full dataset (including both the pre and post-event epochs) and a restricted tail-end subset, beginning at MJD 59800, which was used to assess the asymptotic behaviour of the data while ignoring the complex evolution at early times post-event. To determine the preferred model, we evaluated the difference in BICs:
\begin{equation}\Delta \mathrm{BIC} = \mathrm{BIC}_{\textrm{Exp}} - \mathrm{BIC}_{\textrm{PL}},\end{equation}
where the subscripts denote the BIC scores of the exp. and PL models, respectively. A negative value of
$\Delta \mathrm{BIC}$
indicates that the exponential model is preferred over the power-law model, while a positive value is indicative of the power-law model being favoured.
While neither model perfectly captures the full evolution of the profile change (as seen from the
$\chi^{2}_{\mathrm{r}}$
values in Figure 4), the PL model is predominantly supported in Stokes I, particularly in the tail-end dataset. We summarise the
$\Delta \mathrm{BIC}$
across sub-bands, polarisation, and data sub-samples in Figure 5.
Difference in BIC scores indicating model preference for either the exponential (Exp) or power law (PL) fit to the PC scores for the first principle component of the Stokes I and linear polarisation residuals. A negative
$\Delta\mathrm{BIC}$
value indicates that the exponential model is favoured, whereas a positive value indicates the power-law model is favoured. The full dataset encompasses both pre and post-event epochs. The tail dataset is a restricted subset that commences at MJD 59800 where the fitted curve appears to flatten. In the Stokes Itail dataset, the power-law model is dominant, whereas support for either model is mixed for linear polarisation. Here,
$\Delta\mathrm{BIC}=\mathrm{BIC}_{\mathrm{Exp}}- \mathrm{BIC}_{\mathrm{PL}}$
. For clarity, we indicate
$\Delta \mathrm{BIC}$
for each data sub-set on each cell.
In Figure 4, along with the PC1 scores themselves, we show the best-fit power-law models as a function of time, across frequency and polarisation. In the lower-frequency sub-bands, the scores and fitted curves are very similar for Stokes I and linear polarisation. However, above
$\sim1\,300$
MHz, the separation between PC1 scores for total intensity and linear polarisation grows, with Stokes Iexhibiting higher PC scores than linear polarisation. Notwithstanding this separation, the post-event PC scores grow in magnitude as a function of frequency both for Stokes I and linear polarisation.
This behaviour is reflected in the top two panels of Figure 6, which shows the amplitude and derived half-life from both the exponential decay and power-law model fits to the PC scores as a function of frequency. The model amplitudes show a gradual increase with frequency, though the increase is non-monotonic. Similarly, the derived half-lives decline gradually (but non-monotonically) with frequency. The parameters follow similar trends across the different polarisations and models, although with offsets between them.
Evolution of fitted parameters to Stokes I and linear polarisation PC scores as a function of frequency for the exponential-decay (Exp) and power-law (PL) models (see equations 1 and 2). First Panel: model amplitude A; Second Panel: half-life decay times
$\tau_{1/2}$
; Third Panel:
$\Delta B$
– the difference between the post-event and pre-event baselines; Fourth Panel: S/N of
$\Delta B$
, with the dashed blue line indicating the
$3\sigma$
significance threshold. Shaded regions around each line represent the uncertainties on the fit parameters. Black (light blue) indicates the model fits for the Exp. model in Stokes I (linear polarisation), while red (purple) indicates the model fits for the PL model in Stokes I (linear polarisation).
Both our exp. and PL models include a free parameter
$B_{\mathrm{post}}$
to account for a non-zero asymptotic offset from the pre-event baseline value
$B_{\mathrm{pre}}$
. We measure the absolute value of this offset as
$\Delta B = |B_{\mathrm{pre}} - B_{\mathrm{post}}|$
, and assess the significance of this offset by quantifying the offset S/N as
$|\Delta B|/\sigma_{\Delta B}$
where
$\sigma_{\Delta B}$
is the uncertainty of the offset
$\Delta B$
. We adopt a threshold of 3 as being indicative of a permanent reconfiguration of the profile. We show the behaviour of
$\Delta B$
in the lower two panels of Figure 6. We found that a significant non-zero offset is consistently preferred under the exponential model. However, as discussed above, the power-law model is predominantly supported in Stokes I. In the case of the power-law fits, only some sub-bands showed significant evidence for a non-zero offset: in Stokes I, only sub-band D and sub-band E showed offsets exceeding 3
$\sigma$
significance (
$3.8\sigma$
and
$3.4\sigma$
respectively), whereas in linear polarisation, sub-bands B – E showed offsets with significances ranging from
$3.2$
–
$5.2\sigma$
. In light of these mixed results, it is unclear whether the profile, overall, is fully recovering to its pre-event state, or whether some components may not fully recover. Ongoing long-term monitoring will be required to answer this question more comprehensively.
3.2 Polarimetric analysis of the profile change
Pulsar emission is intrinsically polarised and directional, which results in complex, frequency-dependent polarisation profiles (Xilouris et al. Reference Xilouris1998; Dai et al. Reference Dai2015; Serylak et al. Reference Serylak2021; Oswald et al. Reference Oswald2023b). A key feature of the polarisation in some pulsars is the presence of OPMs, which are thought to be generated as a result of two different modes of wave propagation in the pulsar magnetosphere: the ordinary (O) mode – in which the wave propagation runs parallel to field lines, and the extraordinary (X) mode – in which wave propagation is orthogonally aligned to the orientation of field lines. The intrinsic emission contains components of both these modes. As the waves travel through the birefringent plasma, the O and X mode propagate with different velocities, resulting in a relative phase delay. Upon leaving the plasma, these components recombine and the observed dominant mode depends on the magnitude of this phase delay. As a result, the observed PA can exhibit sudden 90-degree jumps, as one mode becomes dominant over the other across the pulse longitudinal phase. For a comprehensive summary on OPMs, see Gangadhara (Reference Gangadhara1997), von Hoensbroech et al. (Reference von Hoensbroech, Lesch and Kunzl1998), Petrova (Reference Petrova2001), Karastergiou (Reference Karastergiou2009), Noutsos et al. (Reference Noutsos2015), Dyks et al. (Reference Dyks2016), Wang et al. (Reference Wang, Wang and Han2014), Dyks (Reference Dyks2017), Oswald et al. (Reference Oswald, Karastergiou and Johnston2023a).
Like the total-intensity profiles, individual pulses exhibit variability, but the integrated polarisation and PA profiles in MSPs are very stable – with the boundaries of the two modes clearly distinguished. The OPMs also show frequency dependence in their location and extent across the profiles (Smits et al. Reference Smits, Stappers, Edwards, Kuijpers and Ramachandran2006). The superposition and coherent/incoherent mixing of both the O and X mode emission is thought to lead to depolarisation (Stinebring et al. Reference Stinebring, Cordes, Rankin, Weisberg and Boriakoff1984), while the partially coherent mixing of modes has been proposed as an origin of pulsar circular polarisation (Oswald et al. Reference Oswald, Karastergiou and Johnston2023a). Fractional polarisation can therefore provide insights into the overall polarisation profile configuration and the dominance of either of the OPMs across the pulse phase: a high fractional linear polarisation indicates the dominance of one of the modes, whereas low or zero fractional polarisation may indicate mode mixing.
3.2.1 Variability in OPM intensities
J1713 exhibits several OPM transitions, which demarcate four linear polarisation sub-components that we annotate in Figure 1. To analyse how the profile change event affected these components, we first measure the peak normalised flux density of the profile components corresponding to these two OPMs over time across all UWL sub-bands. These are presented in the left column of Figure 7. Specific patterns emerge in the behaviour of the OPM sub-components following the profile change event. The first OPM sub-component (corresponding to negative PA values in Figure 1), for example, shows a significant increase in flux density immediately after the event across all sub-bands, before undergoing an exponential-like decay. The second OPM sub-component (corresponding to positive PA values) is more stable, but exhibits a small but visible amplitude increase immediately after the event. We also measured the ratio of the peak flux densities of the two OPM components over time, for each frequency sub-band. As illustrated in the right column of Figure 7, the profile change event produced a pronounced spike in the flux density ratio, followed by a gradual decay.
Left Column: normalised flux densities of the two profile peaks located in the central phase bins of the linear polarisation profile that represent the two OPM modes, plotted as a function of time and in different frequency bands. The red points indicate the first mode, while blue points indicate the second mode. Right: Flux density ratio of the same two central OPM components. The vertical orange line represents the profile change date (MJD 59320/59321).
3.2.2 Variability in linear polarisation PAs
We examined PA variability using the ephemeris-aligned profile dataset. Figure 8 presents this evolution over time, where we present the PA across OPM components 1.1, 2.1, 1.2 for each sub-band. Notably, the first sub-component (1.1) remains remarkably stable across all epochs, and in particular, its right-hand boundary, where the mode transition to sub-component 2.1 occurs does not appear to shift significantly in pulse phase, nor do the PA values themselves vary. The fact that the PAs do not appear to vary in this leading edge of the profile, while the (linearly polarised) flux density dramatically varies, highlights that this leading edge of the profile is dominated by one orthogonal mode of emission. However, the PA boundary between 2.1 and 1.2 shows significant variation following the profile change event. Specifically, PAs on the left boundary of 1.2 extend leftward, reflecting a contracted phase width for sub-component 2.1. The magnitude of this leftward incursion of component 1.2 into 2.1 increases with frequency, demonstrating a frequency-dependent polarimetric effect of the profile change event.
Time-series of the PA of the central profile sub-components (OPM 1.1 (light red), OPM 2.1 (light blue), OPM 1.2 (dark red)) over the 8 frequency sub-bands, which are labelled in the top left of each sub-figure. We have flagged out PAs in phase bins where the Stokes Isignal-to-noise was lower than 2. The black dashed line indicates the profile change event date. The colourmap spans from
$-70^\circ$
(dark red) to
$70^\circ$
(dark blue).
The impact of this shift in OPMs is also evident in Figure 9, in both the Stokes I and linear polarisation profiles immediately after the event. In Stokes I, profile component C narrows significantly, consistent with the contraction of OPM component 2.1. Profile component D, associated with OPM component 1.2, shows frequency-dependent variability, with the formation of a flat, ‘L-shaped’ plateau-like structure in the lower sub-bands and more gradual gradients in the higher sub-bands. In the linear polarisation profile, the central component narrows only slightly, but the higher sub-bands corresponding to OPM component 1.2 increase in flux density, matching the amplitude seen in lower sub-bands prior to the event, no longer exhibiting the clear separation of frequency-dependent amplitudes for this component.
Top row: PA as a function of pulse phase in the four highest sub-bands, showing template data in black, an early post-event epoch (
$\sim$
50 days after the event; MJD 59368) in red, and a late epoch (
$\sim$
two years post-event; MJD 60049) in orange. A distinct modal bridge emerges between sub-components 2.1 and 1.2, suggesting a complex interaction between the OPMs. Sub-component 1.1 remains stable throughout the profile change event. Bottom row: Ellipticity angle (EA) as a function of pulse phase for the same four sub-bands, showing the template in black, the same early post-event epoch (MJD 59368) in red, and the later post-event epoch in (MJD 60049) in orange.
Figure 9 illustrates the PA sub-components for the four highest sub-bands, comparing the template (black) with the two post-event reference epochs shown in Figure 1: MJD 59368 (red; 210603) and MJD 60049 (orange; 230415), along with the ellipticity angle (EA), defined as
$1/2 \arctan 2(V/L)$
. Examining each sub-component individually confirms that the first component (1.1) remains stable across both time and frequency. In contrast, the third component (1.2) exhibits a pronounced deviation from the template at both epochs, forming a ‘modal bridge’ (i.e. a phase resolved OPM transition) in which the PA values extend leftward and upward toward the second component (2.1), which suggests a complex interaction between the two modes of emission. At similar pulse phases, the EA decreases below
$-22.5^\circ$
, indicating that
$|V| \gt L$
in these regions. The second component itself displays a clear contraction, immediately post-event. These results represent the first evidence of a modal bridge forming after a profile change event in any MSP. This modal restructuring occurs in pulse phase bins where the profile shape itself is altered, indicating a direct link between OPM phenomenology and the overall profile changes. Additionally, these features persist well into the later post-event epoch (MJD 60049), highlighting the long-term nature of this event.
We also note that Figure 1 suggests an apparent shift of the modal bridge towards lower frequencies at later epochs. However, analysis across a wide range of post-event epochs shows that this behaviour is not consistent, with several later epochs exhibiting modal bridge structures in both high and low frequency bands. Our ability to probe this behaviour in more detail is limited by low S/N in many epochs, and we therefore suggest detailed investigation of this effect should be carried out with more sensitive telescopes like MeerKAT or FAST (Miles et al. Reference Miles2025a; Xu et al. Reference Xu2025).
3.2.3 Variability in fractional polarisation
We also investigated the behaviour of fractional polarisation over time, examining phase-averaged OPM profile components over frequency, and the profile overall in a fully phase-resolved fashion. Figure 10 presents the fractional linear polarisation as a function of pulse phase across each sub-band for the template (black), and two high S/N post-event reference epochs: MJD 59368 (red) and MJD 60049 (orange).
Fractional linear polarisation over central bins of pulse phase as a function of both frequency and time. The black curve represents the pre-event template, the red curve represents our first reference epoch (MJD 59368), and the orange curve represents our second reference epoch (MJD 60049). In the bottom left sub-plot, we annotated the sudden increase in
$L/I$
in this region of the profile, as well as the location of the modal bridge.
Fractional linear (squares) and fractional circular polarisation (crosses) per OPM component as a function of both frequency and time. The black data points represent the pre-event template, the red data points represent our first reference epoch (MJD 59368), and the orange data points represent our second reference epoch (MJD 60049).
The template shows a central peak in
$L/I$
corresponding to the main pulse component (profile component ‘C’ in Figure 1), with steep declines to near-zero flux densities on either side (pulse phase 0.49–0.52) corresponding to the transitions between OPM sub-components. In the early post-event epoch (MJD 59368; red curve in Figure 10), the first depolarisation boundary between OPM sub-components 1.1 and 2.1 remains relatively stable across frequency, with a slight deviation at the highest sub-band (sbH; 3 312 MHz). In contrast, the second depolarisation boundary between sub-components 2.1 and 1.2 shifts leftward in phase following the profile event, with the shift becoming more pronounced at higher frequencies. This frequency-dependent shift aligns with the modal bridge observed in the PAs, particularly around 2 320 MHz (sbG), where we annotate this feature as the component 1.2 bridge (see Figure 9). At these phases, the EA also decreases below
$-22.5^\circ$
, indicating that circular polarisation dominates the polarised intensity. By the later reference epoch (MJD 60049; orange curve in Figure 10), the fractional linear polarisation closely resembles the template, indicating a gradual recovery of the polarisation over time.
We also note the emergence of an additional fractional polarisation feature in the highest three sub-bands around pulse phase
$\sim$
0.40–0.44 (Figure 10), where the MJD 59368 profile (red curve) exhibits a significant increase in fractional linear polarisation relative to the template. We annotate this as the component 1.1 profile spike, corresponding to the emergence of a peak in linearly polarised intensity within phases 0.40 – 0.44 in the highest frequency sub-bands, as seen in Figures 1, 2 and 3. We note that this polarisation feature corresponds with the emergence of a persistent new, faint profile feature most apparent in the higher-frequency sub-bands of both Figures 2 and 3, around the location of profile component ‘A’ labelled in Figure 1. Despite the significant change in fractional polarisation in this phase range, the corresponding linear polarisation PA remains unaffected.
We show the average fractional polarisation as a function of frequency, for both linear and circular polarisation (Stokes V), across different OPM sub-components in Figure 11. We compute the average quantities within the respective pulse phase boundaries of each sub-component, which we defined by the continuous extent of the PA clusters seen in e.g. Figure 1. The results are shown for all sub-bands and three key epochs: the pre-event template (black), and our two reference epochs as highlighted previously in Figures 9 and 10, 50 days (red) and
$\sim$
2 yr (orange) post event. Fractional linear polarisation is indicated with square markers and circular polarisation with cross markers. Components 1.1 and 2.1 (representing the leading edge of the profile) show only mild variability across these examples. In contrast, components 1.2 and 2.2 (on the profile trailing edge) exhibit significant variability in the spectral behaviour of fractional polarisation. For component 1.2, the early post-event epoch shows a non-monotonic frequency dependence in fractional linear polarisation, with most sub-bands eventually recovering to near-template values by the late epoch, though some tension remains in the highest frequency bands. Stokes V also displays notable variation in this component. In component 2.2, while Stokes V remains relatively unaffected, fractional linear polarisation significantly varies immediately after the event and had not yet recovered by MJD 60049.
As shown in Figure 11, we also observe the fractional polarisation of the leading sub-components – corresponding to the first half of the polarisation profile – remains relatively stable during and well after the event. In contrast, the trailing sub-components exhibit significant deviations from the pre-event template, even persisting in the later reference epoch. A similar pattern is seen in Stokes V, which is minimally affected in the leading half of the profile but displays greater variability in the trailing region.
4. Discussion
4.1 The origin of the profile change event
Discrete profile changes are rare amongst the MSP population. A comparable case to the J1713 event considered in this paper is the profile change event in PSR J1643
$-$
1224, where a new profile component appeared in observations taken at
$\sim 3\,000$
MHz and
$\sim 800$
MHz (Shannon et al. Reference Shannon2016; Brook et al. Reference Brook2018), and decayed over approximately four months. However, one key distinction from the PSR J1643
$-$
1224 case is that its new component was unpolarised, whereas the variations we report in the J1713 profile change event exhibit clearly distinguished frequency-dependent variations to all Stokes parameters.
Long-lasting profile variations as a result of a discrete profile change were also reported in PSR J0437
$-$
4715 (Goncharov et al. Reference Goncharov2021), with profile residuals persisting for over a year. Similarly, Brook et al. (Reference Brook2018) also documented profile variability in PSR J2145
$-$
0750 and PSR B1937
$+$
21, although these were attributed in part to instrumental systematics, such as polarisation miscalibration, RFI, as well as signal propagation effects such as scintillation and scatter broadening. Outside of these examples, profile change events in the MSP population that exhibit such pronounced variations in the profile shape structure and multi-year recovery times are non-existent.
While the physical mechanism that triggered the April 2021 profile change remains unknown, the differential response between the two OPMs is suggestive of environmental changes in their magnetospheric origin, either close to the emission region (Cheng & Ruderman Reference Cheng and Ruderman1979) or at higher altitudes (Petrova Reference Petrova2001). This is also supported by the frequency-dependence of the fitted parameters to the PCA scores seen in Figure 6, which show that the amplitude of the profile change increases with frequency, while the recovery timescale declines with frequency. This is inconsistent with line-of-sight propagation effects, which therefore makes an intrinsic magnetospheric origin more favourable. However, pinpointing the exact in situ cause of the event remains challenging, owing to the complex and poorly understood nature of MSP magnetospheres (Xilouris et al. Reference Xilouris1998; Stairs, Thorsett, & Camilo Reference Stairs, Thorsett and Camilo1999; Xu et al. Reference Xu2025).
One possibility is that the event was caused by a sudden reconfiguration of the MSP magnetic field, as suggested previously (Goncharov et al. Reference Goncharov2021; Singha et al. Reference Singha2021a; Jennings et al. Reference Jennings2024; Wang et al. Reference Wang2024). Such a reconfiguration would alter the emission geometry and as such, the pulse profile. It may also alter the relative strengths of the OPMs in different pulse phase regions, therefore shifting the pulse phases where depolarisation occurs as each mode transitions in dominance as seen in Figures 8 and 10. The observed profile shape change would then exhibit temporal evolution as the field returned to its prior configuration over time. This picture is consistent with our observation of broad-band variability in the profile and its OPM components, which we note has not been previously reported in connection to this event or other MSP profile changes. Whilst not directly associated with discrete profile changes, Timokhin (Reference Timokhin2010a) and Yuen & Melrose (Reference Yuen and Melrose2017) highlight examples of the variations to pulsar emissions resulting from changes in the magnetospheric configuration and emission geometry.
Additionally, reconfigurations of magnetic field lines are intrinsically linked to a redistribution of plasma density, since active field lines above polar caps – where pair production occurs and emissions are generated – are determined by the global magnetic field structure (Beloborodov Reference Beloborodov2008; Timokhin Reference Timokhin2010b). As such, a redistribution of plasma density could modify the refractive indices for the O and X modes, altering the degree of mode separation, augmenting the mode mixing and leading to observable changes in the relative dominance of either OPM (Wang, Wang, & Han Reference Wang, Wang and Han2015; Noutsos et al. Reference Noutsos2015). This scenario would be consistent with the non-monotonic frequency-dependence observed in the fractional polarisation for certain regions of the profile (marked by the boundaries of the OPM sub-components, see Figure 11 in particular sub-component 1.2), as birefringence effects scale with observing frequency (McKinnon Reference McKinnon1997). Further, the enhanced ellipticity angle at the pulse phases corresponding to new ‘modal bridges’ detailed in Section 3.2.2 is suggestive of partially-coherent mode combination (Oswald et al. Reference Oswald, Karastergiou and Johnston2023a), consistent with changes in the properties of the magnetospheric plasma.
One postulated driver of changes to the magnetospheric plasma is the evaporation of an asteroid as it enters the pulsar magnetosphere, which may cause profile shape changes and alter the pulsar spin properties (e.g. Cordes & Shannon Reference Cordes and Shannon2008; Brook et al. Reference Brook2014). We have not yet carried out a detailed investigated whether the April 2021 profile change was accompanied by variations in the MSP’s rotational period or period derivative before and after the event. This analysis will be presented in an upcoming paper focused on the timing impacts of the event. To date, there have been no reports of changes in the spin properties of J1713 following the event, nor is there any evidence of a large glitch associated with this event.
As discussed in Section 3 and illustrated in Figures 2 and 3, the profile residuals persist through to the latest epoch in our dataset, irrespective of the alignment method used. The absence of similar residuals in the pre-event data under both alignment approaches suggests that this is not a systematic artifact, but a genuine consequence of the April 2021 profile change. This raises a key question: has the pulse profile of J1713 undergone a permanent reconfiguration, or is it still gradually returning to its original configuration? If the former, this would mark only the second such case among MSPs – after the profile reconfiguration in PSR J1643
$-$
1224 reported by Shannon et al. (Reference Shannon2016). However, as outlined in Section 3.1.1 the results of model fitting to the decaying residuals are not conclusive as to whether the profile has permanently reconfigured or is still recovering. Continued monitoring and analysis will be required to resolve this.
Overall, our findings strongly disfavour external or interstellar propagation-related origins for the April 2021 profile change. We find no evidence for intrinsic emission phenomena such as nulling, mode-changing, or glitching, and find no evidence of signal propagation effects including scattering, pulse broadening, or plasma lensing. Pulse-to-pulse variability (jitter) and scintillation are well-characterised for J1713 (Shannon & Cordes Reference Shannon and Cordes2012; Dolch et al. Reference Dolch2014; Bogdanov et al. Reference Bogdanov, Pruszyńska, Lewandowski and Wolszczan2002; Brook et al. Reference Brook2018) and do not account for the long-duration, structured profile changes observed in this event.
4.2 Implications for pulsar timing arrays and future work
J1713 is a high-priority MSP for all global PTA collaborations. However, the April 2021 profile change event introduced large, frequency-dependent timing offsets (Larsen et al. Reference Larsen2024) that led to the exclusion of post-event data from PTA datasets (EPTA Collaboration et al. Reference Collaboration2023b; Xu et al. Reference Xu2023; Reardon et al. Reference Reardon2023; Agazie et al. Reference Agazie2023a; Miles et al. Reference Miles2025a; Rana et al. Reference Rana2025). The exceptional long-term timing stability of J1713 prior to the event means that development of methods to correct for the profile change and enable continued use of this MSP in these high-precision timing campaigns, such as using stable features of the polarisation, should be a priority. While the J1713 event studied in this work is the largest example of its kind for MSPs, these improved techniques to correct for profile changes will prove valuable in the era of highly sensitive observations with FAST and the SKA, which will be capable of detecting smaller profile changes compared to the current catalogue of MSP profile changes.
Rapid-response follow-up observations to profile change events would have scientific benefits beyond potential corrections to timing – in particular, for understanding the phenomenology of pulse shape changes at early times, and subseqeunt insights e.g. into the magnetosphere dynamics driving them (Weltevrede et al. Reference Weltevrede, Johnston and Espinoza2011; Yuen & Melrose Reference Yuen and Melrose2017; Fisher et al. Reference Fisher2024). Additionally, some MSPs exhibit thermal X-ray hotspots near their magnetic poles, caused by returning currents heating the neutron star surface (Zavlin Reference Zavlin2006; Webb et al. Reference Webb2019). If profile changes are indeed driven by magnetospheric field-line reconfigurations, corresponding variability in these X-ray hotspots may also be detectable, which may present an opportunity for X-ray instruments such as XMM-Newton or NICER (Miller et al. Reference Miller2019; Riley et al. Reference Riley2019; Wolff et al. Reference Wolff2021).
Another important avenue for future investigation is the effect of profile change events on single-pulse behaviour, with a particular focus on their polarimetry. For example, during a serendipitously captured glitch in the Vela pulsar, Palfreyman et al. (Reference Palfreyman, Dickey, Hotan, Ellingsen and van Straten2018) observed a dramatic quenching of the Stokes I flux density at the moment of the event, accompanied by quenched linear polarisation in the pulses immediately post-glitch, along with variability in the pulse shape both before and after the glitch. This highlights the potential of single-pulse observations as a diagnostic tool during sudden state changes.
5. Conclusion
We conducted a comprehensive polarimetric analysis of the PSR J1713
$+$
0747 April 2021 profile change event spanning three years of observations with the UWL receiver on Murriyang, providing continuous bandwidth from 704 to 4 032 MHz. Our study was motivated by the insights that such a wide bandwidth, along with high polarimetric fidelity could provide.
Consistent with previous studies, we conclude that the PSR J1713+0747 April 2021 event resulted from a disturbance in the MSP’s magnetosphere, rather than propagation effects through the interstellar medium. In our case, this is supported by the broad bandwidth and long timescale (half-life
$\sim$
160 days) of the profile change. This is further supported with the distinct polarimetric behaviour, such as a stable PA OPM transition in the profile leading-edge seen in Figure 8, whereas the trailing edge of the profile displays shifts in the location of OPM transitions along with greater variability in the fractional polarisation (Figures 10 and 11). Indeed, the relative intensity of the orthogonal modes (Figure 7) shows significant variations after the profile change, connected to complex changes to the profile structure, including features that appear, decay, or sustain over time (Figures 2 and 3). These changes have a substantial, positive dependence on radio frequency (Figure 5). Similarly, while no significant variability is evident in the circularly polarised profiles at lower frequencies, the profile change event has introduced mild variations in Stokes V above
$\sim 2$
GHz (Section 3.1), inconsistent with cold-plasma propagation effects (Cordes, Shannon, & Stinebring Reference Cordes, Shannon and Stinebring2016).
While the magnitude of the profile change has greatly diminished, it is unclear whether the profile is tracking toward full recovery or has permanently reconfigured; our phenomenological models point to either option for different frequency sub-bands. While the results for Stokes I predominantly point towards an ongoing push towards a full recovery, our results for linear polarisation show mixed support for a permanent reconfiguration. Ongoing long-term monitoring will resolve this ambiguity.
The April 2021 PSR J1713
$+$
0747 profile change event, the largest and longest-lasting observed profile change event on any MSP, may provide valuable insights into these complex objects. While stochastic profile variability is common in the canonical pulsar population, discrete profile change events – whilst rare – have been observed in the MSP population. However, as datasets grow and telescope sensitivities improve over time, we may begin to detect profile changes in other MSPs, albeit on smaller scales. Understanding MSP profile changes as we head into the SKA era, their recoveries, and their impacts on ongoing precision timing experiments, such as those conducted with PTAs, is crucial for the continued study of these phenomena and the objects themselves.
Acknowledgement
We thank Simon Johnston and Lucy Oswald for thoughtful discussions and suggestions which improved this work. Murriyang, CSIRO’s Parkes radio telescope, is part of the Australia Telescope National Facility (https://ror.org/05qajvd42) which is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Wiradjuri people as the Traditional Owners of the Observatory site. We acknowledge the Wallumattagal people as the Traditional Owners of the land where this work was carried out. This paper includes archived data obtained through the CSIRO Data Access Portal (http://data.csiro.au). Work at NRL is supported by NASA. MEL is supported by an Australian Research Council (ARC) Discovery Early Career Research Award DE250100508. Parts of this research were conducted by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), project number CE230100016.
Software: PSRCHIVE (van Straten et al. Reference van Straten, Manchester, Johnston and Reynolds2010; van Straten et al. Reference van Straten2011; van Straten, Demorest, & Oslowski Reference van Straten, Demorest and Oslowski2012), TEMPO2 (Edwards, Hobbs, & Manchester Reference Edwards, Hobbs and Manchester2006; Hobbs, Edwards, & Manchester Reference Hobbs, Edwards and Manchester2006)
Profiles, profile residuals, and analysis scripts used in this paper are available here.
Open access
