An Investigation of the state changes of PSR J2021+4026 and the Vela pulsar

We investigate the high energy emission activities of two bright gamma-ray pulsars, PSR~J2021+4026 and Vela. For PSR~J2021+4026, the state changes in the gamma-ray flux and spin-down rate have been observed. We report that the long-term evolution of the gamma-ray flux and timing behavior of PSR~J2021+4026 suggests a new gamma-ray flux recovery at around MJD~58910 and a flux decrease around MJD~59500. During this epoch, the staying time, the gamma-ray flux difference and spin-down rate are smaller than previous epochs in the same state. The waiting timescale of the quasi-periodic state changes is similar to the waiting timescale of the glitch events of the Vela pulsar. For the Vela pulsar, the quench of the radio pulse was in a timescale of $\sim0.2$~s after the 2016 glitch, and the glitch may disturb the structure of the magnetosphere. Nevertheless, we did not find any evidence for a long-term change in the gamma-ray emission properties using years of $Fermi$-LAT data, and therefore, no long-term magnetosphere structural change. We also conduct searching for photons above 100~GeV using 15-year $Fermi$-LAT data, and found none. Our results provide additional information for the relation between the state change of the gamma-ray emission and the glitch event.


INTRODUCTION
Pulsars are highly magnetized and rapidly rotating neutron stars with a stable rotational period.However, there are two types of timing irregularities: glitches and timing noise.A pulsar glitch is defined as a sudden increase in spin frequency.Although the first glitch was discovered over fifty years ago (Radhakrishnan & Manchester 1969), it remains an area of active research and interest.
In addition to timing irregularity, there are various intriguing phenomena observed in pulsar's emission, such as mode changes, nulling, intermittency and pulse-shape variability.It has been inferred that some of the phenomena intimately connect with and arise from alterations in the structure of the pulsar's magnetosphere.For example, Lyne et al. (2010) report that several mode-changing pulsars, which undergo changes in spin-down states, are accompanied by a clear change in radio pulse profiles.Mode changes in pulsars are evident through alterations in the pulse profile, encompassing changes in the relative intensity, rotational phase, and widths of individual profile components.These transitions are consistently observed in radio pulsars, as documented by studies such as Kramer et al. (2006), Wang et al. (2007), and Shaw et al. (2022).
Mode changes likely signify a clear manifestation of current flow redistribution within the magnetosphere, reflecting alterations in the magnetic field, plasma distribution, or other fundamental elements within this region (Timokhin 2010;Huang et al. 2016).
Different timescales and manifestations of mode changes in radio emission and/or the spin-down rate have been observed.For example, PSR J1326-6700 transits between two radio emission states within a timescale of the order of ∼ 100 spin periods (Wang et al. 2007;Wen et al. 2020), while PSR B2035+36 exhibits one abrupt change in spin-down rate on MJD 52950 over 9 years observation, and its radio emission after the mode change is transiting between two states (Kou et al. 2018).On the other hand, several mechanisms for mode changes have been observationally and theoretically suggested.For example, connections between the glitching process of neutron stars and mode changes have been observed through radio observations in several pulsars, including PSRs J1119-6127 (Weltevrede et al. 2011), and B2021+51 (Liu et al. 2021), in which a new component in radio emission appears or the pulse profile varies following glitches.Keith et al. (2013) reported a high correlation between variations in the radio pulse profile and the spin-down rate of PSR B0740-28, which is observed subsequent to the glitch occurring on MJD 55022.Kou et al. (2018) reported a possible glitch event coinciding with a mode change of PSR B2035+36.
As most glitching processes of neutron stars do not show the ra-diative variation, and some mode changes of pulsars occur without any evidence of accompanying glitches, other physical processes have also been suggested to cause mode changes of the pulsars.Precession of neutron stars, which has been suggested as a potential cause for the quasi-periodic changes in the timing residual of PSR B1828-11, exhibits a mode change lasting ∼ 500 days in one cycle (Link & Epstein 2001;Kerr et al. 2016;Jones 2012;Jones et al. 2017).The impact of an asteroid has been also suggested to explain a significant decrease in spin-down rate and variation of the radio pulse profile that occurred in 2005 September for PSR J0738-4042 (Brook et al. 2014), for which the current spin-down rate still remains lower than the value observed prior to 2005 (Zhou et al. 2023;Lower et al. 2023).Shaw et al. (2022) mentioned that the mode change of PSR B2035+36 is similar to that of PSR J0738-4042.Basu et al. (2022) suggested that the evolution of the local magnetic field and the pattern of the sparking discharge in the polar cap acceleration region are causes of the transition between two or more radio emission states.
Variations in emission and mode changes of pulsars have also been observed and investigated with high-energy observations.In particular, since the launch of the Fermi-Large Area Telescope (hereafter Fermi-LAT) in 2008, the known population of γ-ray pulsars has significantly increased.The all-sky monitoring of the Fermi-LAT enables us to investigate the variability and mode change of the pulsars with the γ-ray data.PSR J2021+4026 is a radio-quiet pulsar with a spin-period of P s ∼ 265 ms, and it is known as the first variable γ-ray pulsar.Allafort et al. (2013) reported a sudden decrease in the γ-ray flux by ∼ 14 % and an increase in the spin-down rate by ∆ ḟ / ḟ ∼ 7 %.Subsequent studies (Ng et al. 2016;Zhao et al. 2017) confirmed that the low γ-ray flux state continues for about three years and returns to the previous high γ-ray flux state.Takata et al. (2020) found that PSR J2021+4026 switches between two states with a timescale of several years, namely, a state with high γ-ray flux and low spin-down rate (i.e., HGF/LSD state) and a state with low γ-ray flux and high spin-down rate (i.e., LGF/HSD state).Wang et al. (2018) reported that no significant change in the spectrum and the pulse profile in the X-ray band, which are probably originated from the heated polar cap, has been reported (but see Razzano et al. (2023)).Besides PSR J2021+4026, Ge et al. (2020) reported the model change in the spin-down rate (∼ 0.4%) for PSR J1124-5916, which is a radio-loud γ-ray pulsar, using 12 years Fermi-LAT data.Although PSR J1124-5916 is a glitch pulsar, the correlated glitch events at the mode changes were not confirmed in the Fermi-LAT data.In X-ray bands, Marshall et al. (2015) found a sudden and persistent increase (∼ 36%) in the spin-down rate of the radio-quiet pulsar PSR B0540-69, utilizing data from RXT E and S wi f t telescopes.
With the timing solution obtained from radio observation or Fermi-LAT data, about 50 Fermi-LAT pulsars show glitch activity, and its long-term timing solutions obtained from Fermi-LAT data have been utilized to constrain the glitch mechanisms (e.g.Gügercinoglu et al. 2022).Lin et al. (2021) investigated the variability of the γ-ray emission from the glitching pulsar, PSR J1420-6048, and found the spectral variation of the γ-ray emission between each glitch.Several γ-ray pulsars with a characteristic spindown age similar to that of PSR J2021+4026 (τ s ∼ 7 × 10 4 years) have displayed glitching activities 1 .Hence, It is not unexpected that PSR J2021+4026 has also experienced the glitch activities, though such an event has not been confirmed from the long-term timing 1 https://www.jb.man.ac.uk/pulsar/glitches/gTable.htmlephemeris derived from the Fermi-data.Since the mechanism of the mode change for PSR J2021+4026 has not been well understood, further investigation on PSR J2021+4026 and studies on other glitching γ-ray pulsars will be desired to examine the influence of the glitch activities on the γ-ray emission properties.
The Vela pulsar (PSR B0833-45/PSR J0835-4510) is associated with the Vela supernova remnant, and was discovered over 50 years ago in 1968 (Large et al. 1968).The Vela pulsar is one of the brightest gamma-ray sources in the sky, with a spin period ∼ 89 ms and spin-down rate ∼ 1.25 × 10 −13 ss −1 .It is known as the first pulsar that was observed to undergo glitch (Radhakrishnan & Manchester 1969) and 24 events have been reported to-date 1 .Among all glitching pulsars, the Vela pulsar shows comparatively frequent glitching activity with a waiting time of ∼3 years (Dodson et al. 2007), and a relatively large glitch size of △ν/ν ∼ 10 −6 (Espinoza et al. 2011).Interestingly, Palfreyman et al. (2018) reported a sudden change in the radio pulse shape coincident with the glitch activity on 2016 December 12; an evolution of the pulse profile and polarization was observed in several rotation cycles of the Vela pulsar.The results indicate that the glitch event of the Vela pulsar has an impact on the structure of the magnetosphere.Since the Vela pulsar is one of the brightest gammaray pulsars, it is worthwhile to investigate whether the glitch activity causes any variation or state change in the gamma-ray emission.
In this paper, we study the temporal evolution of the gammaray emission properties for two bright gamma-ray pulsars, PSR J2021+4026 and Vela pulsar.We created a timing solution for each state of PSR J2021+4026 and for each inter-glitch interval of the Vela pulsar.We compared the pulse profiles and spectra in different time intervals.Section 2 describes the data reduction process for the two pulsars observed by Fermi-LAT.In section 3, we present the results of our data analysis and report a new state change of PSR J2021+4026.In Section 4, we provide a discussion on the implications of our results.

Fermi-LAT data
Fermi-LAT is a gamma-ray imaging instrument that scans the whole sky approximately every three hours, covering the energy band from ∼ 20 MeV to 300 GeV (Atwood et al. 2009).We selected Pass 8 data in the energy band of 0.1-300 GeV.We choose the position to be R.A.=20 h 21 m 30 s .48,decl.=40o 26 ′ 53 ′′ .5 for PSR J2021+4026, and R.A.=08 h 34 m 00 s .00,decl.=45o 49 ′ 48 ′′ .0 for the Vela pulsar.To avoid Earth's limb contamination, we only included events with zenith angles below 90 degrees.We limited our analysis to events from the point source or Galactic diffuse class (event class = 128) and used data from both the front and back sections of the tracker (evttype = 3).For our spectral analysis, we constructed a background emission model that incorporated both the Galactic diffuse emission (gll_iem_v07) and the isotropic diffuse emission (iso_P8R3_SOURCE_V3_v1) provided by the Fermi Science Support Center.For each bin of spectra and flux evolution, we refit the data by the binned likelihood analysis (gtlike) and estimate the flux.
To study the evolution of the gamma-ray flux above 0.1 GeV, we divided the data into 60-day time bins.The contribution of each background source is calculated with the spectral parameters obtained from the entire data in each bin.Then, we refit the data by the binned likelihood analysis (gtlike) and estimate the flux.
We also performed the likelihood analysis of the Vela pulsar in the time range of 2008 August to 2022 December.We use the power-law with an exponential cutoff model as mentioned above, and the data are divided into seven time segments separated by glitches: MJD .By creating timing ephemeris for each interglitch interval, we study the evolution of the gamma-ray flux and examine the pulse profile in different energy bands.This will allow us to determine whether there are any significant changes in the pulsar's high-energy emission that is associated with glitch events.

Timing analysis
To ensure that we include the majority of significant source photons, we consider the photon events within a 1 • aperture centered at the targets.We assign phases to these photons using the Fermi plug-in for TEMPO2 (Hobbs et al. 2006;Edwards et al. 2006).We barycentred the photon arrival times to TDB (Barycentric Dynamical Time) using the JPL DE405 Earth ephemeris (Standish et al. 1998).This was done using the "gtbary" task, which provides accurate measurements of the positions and velocities of Solar System.To obtain the timing ephemeris, we use the Gaussian kernel density estimation method (de Jager et al. 1986) provided in Ray et al. (2011) to build an initial template.We then use this template to cross-correlate with the unbinned geocentered data to determine the pulse time-ofarrival (TOA) for each pulse.Once the pulse TOAs are obtained, we can fit them to an original timing model assumed from the semiblind search, the timing results of PSR J2021+4026 and Vela pulsar are shown in Appendix A.
To obtain the temporal evolution of the spin-down rate, we divide the data into several ten-day time bin and determine the first derivative of the frequency for each time bin.We also create the long-term ephemeris for the two states of PSR J2021+4026 and for each interglitch interval between two glitch events for the Vela pulsar.

Temporal evolution
The top and bottom panels in Figure 1 show the flux temporal evolution and first-time derivative of the spin frequency ḟ , respectively; each data point is obtained with the 60-day Fermi-LAT observation.In addition to the previous three events of the state changes reported by Allafort et al. (2013) and Takata et al. (2020) (vertical black dashed lines in Figure 1), we find new state change from LGF/HSD to HGF/LSD around MJD 58910.Prokhorov & Moraghan (2023) also reported the transition to a high-flux state in June 2020 but they did not discuss any changes in the spin-down rate associated with this event and from HGF/LSD to LGF/HSD around MJD 59510.Since we cannot determine a precise onset time of the state change, we indicate an approximate location with the red vertical dashed line in Figure 1.Hereafter, we denote the newly-found states as LGF/HSD 3 and HGF/LSD 3, as indicated in Figure 1.
As Figure 1 indicates, the time durations of LGF/HSD 2 and HGF/LSD 3 are shorter than those in the previous states.For example, LGF/HSD 2 continues about 840 days, which is significantly shorter than about 1140 days (> 1140 days) for HGF/LSD 2 (HGF/LSD 1).We also find that the average flux of the new state is different from the previous values.The flux ∼ 1.12(1) × 10 −6 photon cm −2 s −1 of new HGF/LSD 3 is slightly smaller than 1.15−1.16×10−6 photon cm −2 s −1 of the previous HGF/LSD 1 and 2. The flux level of new LGF/HSD 3 is higher than that in the previous two states of LGF/HSD 1 and 2. The flux changes from HGF/LSD to LGF/HSD in previous two events were ∼ 12 − 14% (ref.Table 1), while the fractional flux change of the new event is only ∼ 4%.
A sudden change of the first-time derivative of the frequency, ḟ , with a time scale of < 10 days was observed at the previous state changes from HGF/LSD to LGF/HSD (Allafort et al. 2013;Takata et al. 2020).For new event from HGF/LSD 3 to LGF/HSD 3, on the other hand, a sudden change in timing behavior could not be confirmed, and the ḟ shows a more gradual evolution to HSD as indicated in the bottom panel of Figure 1.We obtained the spin-down rate at each point with the data of a 60-day time bin.Neglecting high-order timing noise, we describe the evolution of timing solution using a linear expression as f (t) = f (t 0 ) + ḟ (t 0 )(t − t 0 ), where, f and ḟ are spin frequency and frequency derivative, respectively.The middle point of each 60-day time bin is chosen to determine the reference time (t 0 ), the uncertainty is determined from the Fourier width, 1/70 days, and the best frequency and spin-down rate in each bin are determined by the results to gain the most significant detections using H-test (de Jager & Büsching 2010).In order to examine the track of the long-term behavior, we calculate the average time derivatives for each state (i.e., horizontal dashed lines in the bottom panel of Figure 1).We find that the average value of ḟ = −7.81(5)×10−13 Hz s −1 in new HGF/LSD 3 is smaller than those in previous HGF/LSD states, while ḟ = −8.13(5)× 10 −13 Hz s −1 in new LGF/HSD 3 is consistent with the previous LGF/HSD states within errors.Figure 2 summarizes the evolution of the spin frequency from HGF/LSD 3 to LGF/HSD 3. The spin-down rate at each point is calculated with the data of a 70-day time bin, and the starting time difference between two neighbor points is 10 days; namely, two neighbor points overlap the 60-day data.In this figure, we can see a dramatic change of the spin-down rate at around MJD 59510-59600.With the current uncertainty of the data, on the other hand, we cannot confirm the existence of a glitch.If there were a glitch, the frequency jump would be smaller than △ f < 10 −7 Hz, which is consistent with the previous events from HGF/LSD to LGF/HSD (Allafort et al. 2013;Takata et al. 2020).

Timescale of state change
In this section, we investigate the correlation between the spin-down rate and the gamma-ray flux of PSR J2021+4026 from 2008 to 2023 by utilizing the discrete correlation function (DCF).The DCF mea- sures correlation functions without interpolating in the temporal domain (Edelson & Krolik 1988).The unbinned discrete correlation function is defined as: where a i and b j represent the data lists of gamma-ray flux and spindown rate, respectively, in each time bin of Figure 1, and a and b represent the average values of the a i and b j , respectively.In addition, σ a and σ b are the standard errors of a and b, respectively.We calculate the time-lag of each pair (a i , b j ) and create the group of the pairs for every 10 7 seconds of the time-lag (0 < △t 1 < 10 7 s, 10 7 s < △t 2 < 2 × 10 7 s, ...).For each bin of the time-lag, we obtain number of the pairs, M, and calculate the average of UDCF, where a standard error for the DCF is defined as, In Figure 3, we present the DCF curve, which exhibits the maximum correlation coefficient of ∼ 0.32 at a time lag of zero seconds.This indicates that the derivative of frequency changes its state simultaneously with the LAT flux state.The values of the DCF show periodicity with a period of 2 × 10 8 seconds, which is approximately 6.5 years.This result is consistent with that of Takata et al. (2020), which suggests that the time interval between the starting of a specific HGF/LSD is about 6-7 years.

Phase-averaged spectrum and pulse profile
In Figure 4, we present the observed spectra for the new HGF/LSD 3 (black dots) and LGF/HSD 3 (red dots) states.The solid lines show the fitting function using Equation 1. Table 2 summarizes the parameters of the best fitting function.As indicated by the timing evolution (Figure 1 and Table 1), the flux change from HGF/LSD 3 to LGF/HSD 3 is only ∼ 4%, which is much smaller than 12 − 14% of the previous two events.The previous state changes from HGF/LSD to LGF/HSD make the spectrum softer; an increase of γ 1 and a decrease of the cut-off energy were observed.From the HGF/LSD 3 to LGF/HSD 3, both the power-law index γ 1 and the cut-off energy increase.
To investigate the changes in the pulse profile after each state change, we accumulated photons of each new state to derive a global 1 Average flux in each state (10 −6 photon cm −2 s −1 ).
2 Difference between the average flux of each state and averaged flux of whole data (10 −6 photon cm −2 s −1 ). Figure 2. Difference between the observed frequency and predicted frequency from the ephemeris ( ḟ , in the last second column (HGF/LSD 3) in Table A1) around MJD 59510.The observed frequency is determined with the data of the 70 day time bin, and two neighbor points overlap the 60-day data.The uncertainty is determined from the Fourier width, 1/70 days. 1 Energy flux of each state in units of 10 −10 erg cm −2 s −1 .ephemeris for HGF/LSD 3 and LGF/HSD 3 (Table A1), and then folded the corresponding pulse profiles (as shown in Figure 5).We fit the obtained pulse profiles with two Gaussian functions.Table 1 summarizes the parameters of the fitting functions for each state.In the previous state changes from HSG/LSD state to LGF/HSD, the ratio of height for the first pulse (small pulse) to the height of the second pulse (main pulse) decreased from about 50% to 25%, as Table 1 indicates.From HGF/LSD 3 to LGF/HSD 3, we can see a decrease in the ratio of the pulse height, but the magnitude of the change is smaller that in the previous cases.This would be consistent with the smaller flux change in the new event compared to the previous cases.LGF/HSD 3. The solid lines show the best fit models for each spectrum according to the 'PLSuperExpCutoff2' model (see Equation 1).

Long-term behaviors
After the launch of the Fermi telescope, six glitches of the Vela pulsar have been observed.We search for the semi-permanent state changes of the GeV emission triggered by glitches.Figure 6 presents the evolution of the spin-down rate and the gamma-ray flux of the Vela pulsar.We divide the entire Fermi-LAT data into seven time intervals, which are bounded by the time of the glitches.The top panel of Figure 6 shows the temporal evolution of the gamma-ray flux, and the bottom panel shows the evolution of the spin-down rate.Timing ephemeris and the emission properties for each inter-glitch interval are summarized in Table A2 and Table 3, respectively.Despite the frequent glitches of the Vela pulsar, we do not confirm any significant state change in the spectral properties associated with the glitch.In the bottom panel of Figure 6, the flux shows a small fluctuation with an amplitude of several percentages.However, such a fluctuation can be explained by the random fluctuation observed in the data.Table 3 summarizes the information of the time-averaged spectrum for each inter-glitch interval.We can find that the spectral properties of the energy flux, the spectral index (γ 1 ) and the cut-off energy in the different time intervals are consistent with each other within a limited percentage of errors.
In Figure 6, we determine the average flux with the 1 σ uncertainty for each inter-glitch interval (labeled by the black dashed line and light-blue region).Since the flux uncertainties of each interglitch interval are encompassed within the range of the long-term average flux, we conclude that there is no significant variation of the flux in different inter-glitch intervals.We also note that one or two data points in each inter-glitch time interval may exhibit notably high or low fluxes, and it is likely caused by random fluctuations.As indicated by Table 3, no significant evolution of spectral properties in different inter-glitch intervals is found.

Pulse profile at different time intervals
The brightness of gamma-ray emission enables us to investigate the temporal evolution of the structure of the pulse profile.In order to explore the long-term evolution trend of the pulsed shape, we accumulate ∼ 10 4 photons in the energy range of 0.1-300 GeV to fold a series of pulse profiles using the ephemeris of each inter-glitch inter-val, as given in Table A2.We employ a fitting approach on the pulse profile using a combination of four Gaussian functions, which can better fit a plateau between two main pulses, as illustrated in Figure 7.The correlation coefficient of 0.99, derived from the optimal fit employing four Gaussians, surpasses the correlation coefficient of 0.86 obtained when using only three Gaussians to fit the data.With a series of the pulse profiles generated with 10 4 photons, we then determine the temporal evolution of the phase separation between pulse 1 and pulse 2, the Gaussian width of pulse 1 and the Gaussian width of pulse 2 as shown in Figure 8.We find that there is no significant semi-permanent change shown in the long-term evolution of the pulse profile.The glitch of the Vela may still cause some disturbances of the magnetosphere, as indicated in the radio data (Palfreyman et al. 2018).However, its effect may be on shortterm timescale, or is so small that it is not seen in the gamma-ray emission properties averaged over a time scale much longer than 10 days.

Searching for >100
GeV photon with Fermi-LAT data Leung et al. (2014) searched for the pulsed emission above 50 GeV using about 5.2 years of Fermi-LAT data and reported 5 photons that probably originated from the Vela pulsar.They also report the highest energy photon about 208 GeV with a significant level of 2.2 σ.
The results of H.E.S.S.-II (H.E. S. S. Collaboration et al. 2018) also suggest that the spectrum of the pulsed emission from the Vela pulsar extends beyond 100 GeV.
We revisit the search for very high-energy gamma-ray photons with Fermi-LAT data using updated Fermi-LAT gamma-ray catalog (PASS 8 data), galactic diffuse emission (gll_iem_v07.fits)and isotropic background emission (P8R3_SOURCE_V3) models, while Leung et al. (2014) used the reprocessed Pass 7 "Source" class (P7REP_SOURCE_V15 IRFs).Following Leung et al. (2014), firstly, we divide the Fermi-LAT data obtained within MJD 54710-60000 into 50-300 GeV with a ROI (region of interest) of 4-degree in radius.To select photons originating from the pulsar, we used the gtsrcprob tool to calculate the probability of each photon within ROI associated with the pulsar.We also calculated the probabilities of the photons coming from the Vela X (P PWN ) and the Galactic diffuse emission (P GAL ).We selected only photons with P Vela > max(P PWN , P GAL ).Subsequently, the same methodology was applied to each state, utilizing the respective ephemeris obtained from the timing analysis.
Figure 9 shows the light curve in different energy bands from panel (a) to (c), the panel (d) shows the weighted light curve above 50 GeV.The results in the MJD 54686-56583 in comparison to Leung et al. (2014) are depicted in panel (e) of Figure 9, which illustrates the energy distribution of photons.We detect 7 photons with an energy larger than 50 GeV shown in Table 4.The arrival times of two photons, observed at MJD 56149 and MJD 56317, are consistent with the results reported by Leung et al. (2014).However, the energy of the photon detected at MJD 56317, ∼ 54.6 GeV, is significantly smaller than ∼ 79.5 GeV of the previous result.Although we employ the same method as Leung et al. (2014) to investigate the emission in different time ranges, we do not confirm any photons above 100 GeV.

DISCUSSION AND SUMMARY
We have investigated the temporal evolution of the characteristics of the GeV emission from two young pulsars, PSR J2021+4026 and  av Average flux in each state (10 −6 photon cm −2 s −1 ).
b Ratio of amplitude.the Vela pulsar.For PSR J2021+4026, we confirm new state changes from LGF/HSD to HGF/LSD and from HGF/LSD to LGF/HSD, but the magnitude of the change in the gamma-ray flux and the spindown rate is smaller than those seen in previous state changes.The pulse shape change is also small compared to previous events.Concerning the state change of the emission and spin-down properties, we speculate that the change of global electric current circulating the magnetosphere could lead to corresponding state change of the emis-sion and spin-down properties, and a larger (smaller) electric current produces larger (smaller) spin-down rate in LGF/HSD (HGF/LSD) state.We expect that the change of gamma-ray flux is due to the change of the size of the particle acceleration/emission regions, and HGF/LSD (LGF/HSD) has a larger (smaller) size of the gap.There is a tendency in the previous states that the HGF/LSD state has a larger cut-off energy in the spectrum and the width of the pulse profile is wider (Takata et al. 2020).This tendency can be understood if the HGF/LSD has a larger acceleration/emission region.Although it is expected that the magnitude of the electric current and the size of the acceleration/emission region are related (Hirotani 2006), a positive or negative correlation will depend on the location of the acceleration/emission region and the geometry of the magnetosphere (Takata et al. 2006).For the most recent state change from HGF/LSD to LGF/HSD, the small change in the spin-down rate implies the change of the electric current is smaller, which results in a small change in the structure of the acceleration/emission region.Takata et al. (2020) proposed that the crust cracking process triggers the mode change of PSR J2021+4026, and it causes a sudden change in the magnetic field structure of the polar cap region.According to this model, the timescale between each mode change from HGF/LSD state to LGF/HSD state corresponds to the timescale to accumulate the magnetic stress inside the crust, and the waiting timescale in LGF/HSD state is the relaxation timescale during which the polar cap structure returns to its original state.Jones (2012) suggested a possible correlation between the neutron star's spin period and the waiting timescale of mode change of the transient radio pulsar, and argued the idea of free precession of neutron stars with strained crusts.In this model, the capacity of the acceleration of particles within the magnetosphere was linked to different phases of precession.The precession model has been widely applied to the observed periodic variations of the spin-down rate of PSR B1828-11 (Link & Epstein 2001;Ashton et al. 2017).For PSR J2021+4026, we found that the waiting time in LGF/HSD 2 (∼ 840 days) and HGF/LSDF 3 (∼ 530 days) were shorter than the previous corresponding states (∼ 1140 days), indicating that the mode change is quasi-periodic event.Such a quasi-periodicity may be incompatible with the free-precession model, for which the stable periodic variation of the spin-down rate will be expected.
Among the transitioning radio pulsars, PSRs J0738-4042 and B2035+36 have shown a mode change with relatively large variations in the spin-down rate, |∆ ḟ |/ f ∼ 7% and 12%, respectively.The magnitude of the spin-down rate change of the two pulsars is similar to |∆ ḟ |/ f ∼ 7% of PSR J2021+4026.Moreover, the waiting time to the next mode change of three pulsars will be of the order of or longer than ∼ 10 years; it is observed once or twice mode changes of PSR J0738-4042 over 50 years (Lower et al. 2023) and once of PSR B2035+36 over 9 years (Kou et al. 2018).Brook et al. (2014) argued the mode change of PSR J0738-4042 with an asteroid encounter model (Cordes & Shannon 2008).In this model, the aster-oid with neutral charge can approach to the pulsar magnetosphere, and it is evaporated and ionized by an X-ray irradiation of the pulsar.The ionized gas that migrates within the magnetosphere changes the torque exerted on the neutron star and the emission properties.Due to the size of the frequency jump and waiting time in several years, the asteroid encounter may be an attractive model as the origin of the mode change of PSR J2021+4026.
Since about 50 glitching Fermi-LAT pulsars have been confirmed 1 , it is not unexpected that PSR J2021+4026 has experienced the glitch process of the neutron star.In this paper, we focus on the Vela pulsar, which is both a glitching and bright γ-ray pulsar, we analyzed Fermi-LAT data of the Vela pulsar to investigate the influence of glitch on its γ-ray emission.Palfreyman et al. (2018) reported sudden changes in the radio pulse shape coincident with the 2016 glitch event of Vela pulsar, indicating the glitch affected the structure of the magnetosphere.Bransgrove et al. (2020) suggested that the glitch launched the Alfven wave, which subsequently caused the high-energy radiation and the electron-positron pair creation.They suggest that the created pairs quench the radio emission at the 2016 event of the Vela pulsar.Although the enhancement of the gamma-ray flux is expected during propagation of the Alfven wave in the magnetosphere, the predicted timescale of the existence of the wave, ∼ 0.2 s, is too short to be investigated with Fermi-LAT data.We also did not find any evidence for the change of the magnetosphere structure for the Vela pulsar in a timescale of years, as discussed in section 3.2.Consequently, we cannot conclude that the glitch is the trigger of PSR J2021+4026 based on the current study. .From (a) to (c): the folded light curve in the energy range of 0.1-1 GeV, 1-10 GeV, and 10-100 GeV using the photons collected within 1 • radius centered at the source.From (d) to (e): the weighted light curve in the energy of 50-300 GeV and the energy distribution of the photons above 50 GeV, which are collected within 4 • radius and with a source probability of P Vela > max(P PWN , P GAL ) in the same way as Leung et al. (2014) .The data collected during MJD 54686-56583 as an example to find the energy of photons.
In summary, we have examined the evolution of gammaray flux and spin-down rate in two bright gamma-ray pulsars, PSR J2021+4026 and the Vela pulsar.We confirmed new states change from LGF/HSD to HGF/LSD and from HGF/LSD to LGF/HSD observed for PSR J2021+4026, occurring around MJD 58910 and MJD 59510, respectively.We found that the waiting time, the change in the flux and the first time derivative of frequency are smaller than previous events.The change in the shape of the pulse profile is also smaller than those of the previous events.Our results confirm that the state change of PSR J2021+4026 is not regularly repeated but quasi-periodically repeated with the different magnitude of change of the gamma-ray flux and the spin-down rate.We also carried out the timing and spectral analysis of the Vela pulsar to investigate the effect of the glitch on the observed gamma-ray emission properties, since the radio observation indicates that the glitch disturbs the structure of the magnetosphere.However, we did not confirm any significant state change of the gamma-ray emission triggered by the glitch of the Vela pulsar.Finally, using the 15-year Fermi-LAT data of the Vela pulsar to search photons above 100 GeV, we did not find any photons above 100 GeV from the Fermi-LAT data.We still did not find any evidence for the change of the magnetosphere structure in a timescale of years in the 2016 glitch, while the created pairs may quench the radio emission.
Table A2.The ephemerides of Vela before and after each glitch.

Figure 1 .
Figure 1.Top: gamma-ray flux (> 0.1 GeV) evolution of PSR J2021+4026 from 2008 to 2023.The light-blue/red region labels the uncertainty of the average flux in different states in 1 σ confidence interval.Bottom: The spin-down rate of PSR J2021+4026 from 2008 to 2023.Each point was obtained from the data with a 60-day time bin.The black vertical lines represent the state changes reported in Allafort et al. (2013) and Takata et al. (2020).The vertical red dashed lines indicate new state change reported in this study, and the right-most one is determined by the start of the flux drop.The horizontal dashed lines indicate the averaged spin-down rate for each state.

Figure 3 .
Figure 3.The correlation of Fermi-LAT flux and spin down rate of PSR J2021+4026.The data is from Fermi-LAT in the time range of 2008 to 2023.

Figure 4 .
Figure 4.The comparison of spectra at the state of HGF/LSD 3 andLGF/HSD 3. The solid lines show the best fit models for each spectrum according to the 'PLSuperExpCutoff2' model (see Equation1).

Figure 5 .Figure 6 .
Figure5.The comparison of the pulse profile for PSR J2021+4026 at the time of HGF/LSD 1, LGF/HSD 1, HGF/LSD 2, LGF/HSD 2, HGF/LSD 3 and LGF/HSD 3. The pulse profile is generated with photon energy >0.1 GeV.The vertical dotted lines show the phase interval of pulse 1 and pulse 2.

Figure 7 .
Figure 7.The pulse profile of Vela pulsar (MJD 57740-58510) fitted with 3 (left) and 4 (right) Gaussian functions.The black lines show the pulse profile with the weighted photons, while the red lines depict the results of the Gaussian function fitting.The vertical-dashed-balck lines show the phase intervals of pulse 1, pulse 2, and plateau.

Figure 8 .
Figure8.The parameters to describe the pulsed structure in different time intervals of the Vela pulsar.Each data point was determined using at least 10000 photons.The upper panel shows the phase separation between pulse1 and pulse 2 (∆Φ 21 ), respectively, as shown in Figure7.The middle and the lower panels show the Gaussian width of pulse 1 and pulse 2 (W1 and W2), respectively.The units of ∆Φ 21 , W1 and W2 are the phase in one spin cycle.
Figure9.From (a) to (c): the folded light curve in the energy range of 0.1-1 GeV, 1-10 GeV, and 10-100 GeV using the photons collected within 1 • radius centered at the source.From (d) to (e): the weighted light curve in the energy of 50-300 GeV and the energy distribution of the photons above 50 GeV, which are collected within 4 • radius and with a source probability of P Vela > max(P PWN , P GAL ) in the same way asLeung et al. (2014) .The data collected during MJD 54686-56583 as an example to find the energy of photons.
numbers in parentheses denote 1σ errors in the last digit.

Table 1 .
Information of the average flux and the parameters of the pulsed structure for PSR J2021+4026.

Table 3 .
Parameters of spectra of the Vela pulsar at different inter-glitch intervals.

Table 4 .
Energy, arrival time and source probability of >50 GeV photons with P PS R large than max(P PWN , P PGAL ) for the Vela pulsar.