Repeated State Change of Variable Gamma-Ray Pulsar PSR J2021+4026

PSR J2021+4026 (3FGL J2021.5+4026) is one of the brightest radio-quiet gamma-ray pulsars (Abdo et al. 2009) with a spin period of P_s ∼ 0.265 s and a spindown power of L_sd ∼ 3.4 × 10³⁶ erg s⁻¹. This pulsar is also known as the first variable gamma-ray pulsar (Allafort et al. 2013, hereafter Al13). Fermi-Large Area Telescope (Fermi-LAT) observed a sudden gamma-ray flux drop near MJD 55,850 (2011 October) and a change of the observed flux (>100 MeV) from (8.33 ± 0.08) × 10⁻¹⁰ erg cm⁻² s⁻¹ to (6.83 ± 0.13) × 10⁻¹⁰ erg cm⁻² s⁻¹ over a timescale of less than a week. The timing analysis also revealed that the spindown rate changed from −(7.8 ± 0.1) × 10⁻¹³ Hz s⁻¹ to −(8.1 ± 0.1) × 10⁻¹³ Hz s⁻¹. Subsequent studies (Ng et al. 2016; Zhao et al. 2017) found that the gamma-ray flux and the spindown rate did not show a gradual recovery toward the values before 2011 October, and the pulsar stayed at low gamma-ray flux and high spindown rate state (hereafter LGF/HSD state) for about three years after the state change. These flux and timing behaviors suggest that the state change of PSR J2021+4026 is owing to the change of the structure of the global magnetosphere. Around 2014 December, PSR J2021+4026 returned to the high gamma-ray flux/low spindown rate state (hereafter HGF/LSD state) over a timescale of a few months.

PSR J2021+4026 stayed at the HGF/LSD state for about three years since 2014 December. The aperture photometry light curve (time resolution of 30 days) of PSR J2021+4026 provided by the Fermi science team⁹ clearly shows that a gamma-ray flux change occurred around MJD 58,150 (2018 February 1) and that PSR J2021+4026 entered new low gamma-ray flux state. Considering that the previous flux change accompanied a change in the spindown rate, it is likely that there was another state change in the spindown rate near MJD 58,150. In the previous study (Zhao et al. 2017), we investigated the evolution of the gamma-ray flux and spindown rate until MJD 57,700 (2016 November 8). In this study, therefore, we perform timing and spectral analyses for gamma-ray emission of PSR J2021+4026 with the data after MJD 57,700 in Sections 2 and 3. We will discuss the possible mechanism of the state switching of the PSR J2021+4026 in Section 4.

We use the Fermi-LAT Pass 8 (P8R2) data in the energy range of 0.1–300 GeV for the data analysis. To obtain the temporal evolution of the flux and spindown rate, we select the data from 2016 November 8 to 2019 March 1 in a 20° radius region of interest (ROI) centered at the 3FGL J2021.5+4026 position. The light curve of PSR J2021+4026 is obtained at E > 0.1 GeV using the binned likelihood analysis in the Fermi Science Tools. We select events in the class for the point source or Galactic diffuse analysis (event class = 128) and consider the photons collected in the front and back sections of the tracker (evttype = 3). We exclude events with zenith angles larger than 90° to reduce the Earth's albedo gamma-ray contamination.

For the spectral analysis, we input spectral models, which are provided by the Fermi Science Support Center¹⁰ , for the galactic diffuse emission (gll_iem_v06), the isotropic diffuse emission (iso_P8R2_SOURCE_V6_v06) and all the 3FGL catalog sources within 20° of PSR J2021+4026 to account for the spectral contribution. We also include the extended sources, i.e., the Cygnus-X cocoon and Cygnus Loop, near the PSR J2021+4026. We first fit our data over the entire time range and obtain spectral parameters for the input sources. In this study, the spectrum of PSR J2021+4026 is modeled with a power law plus exponential cutoff function of the form

$$\begin{eqnarray}&&{dN}/{dE}\propto {E}^{-{\rm{\Gamma }}}\exp [-(E/{E}_{c})],\end{eqnarray} \tag{ 1 }$$

where E is the energy of the photon, Γ is the photon index, and E_c is the cutoff energy. To obtain the best-fitting spectral parameters, we perform a binned likelihood analysis using the Fermi Science Tools gtlike. We divide the entire energy range (0.1–300 GeV) into 12 energy bins with a logarithmic scale, and best-fit spectral parameters are obtained with the elimination of insignificant sources (<3σ). The 1σ uncertainty of the fitting parameters is provided by gtlike tool, and it is calculated by the method discussed in Mattox et al. (1996).

We divide the entire time range (from 2016 November 8 to 2019 March 1) into 30 day time bins and obtain the long temporal evolution of the flux in E > 0.1 GeV bands. For each bin, the contribution of each background source is calculated with the spectral parameters that are obtained from the entire data. For PSR J2021+4026, we refit the data by the binned likelihood analysis (gtlike) and estimate the flux.

To obtain ephemeris and the pulse profile of PSR J2021+4026 after 2018 state change, we extract the source events within 1° radius centered at the target. The time of arrival of all photons is corrected to TDB (barycentric dynamical time) with the task gtbary supported by Fermi Science Tools.

3.1. Evolutions of the Gamma-Ray Flux and Spindown Rate

Top panel of Figure 1 shows the long-term evolution for the >0.1 GeV flux with 30 day time bins that was used in Al13. In addition to the previous event near MJD 55,850 reported in Al13, an obvious flux change occurred at around MJD 58,100–58,200. The photon flux dropped by about 13 ± 1% from F = (1.29 ± 0.01) × 10⁻⁶ cts cm⁻² s⁻¹ averaged over MJD 57,105–58,120 to F = (1.12 ± 0.01) × 10⁻⁶ cts cm⁻² s⁻¹ averaged over MJD 58,265–58,500. The size of the flux drop for the new event is consistent with 11 ± 1% of the previous event within the error (Table 3).

Zoom In Zoom Out Reset image size

In the previous flux change, a sudden change of the spindown rate was observed. To examine the evolution of the spindown rate at around the epoch of the new flux drop, we divide the Fermi-LAT data after ∼MJD 57,700 into a 70 day time bin, which can provide enough photons to detect the pulsation (Zhao et al. 2017). During the analysis, we notice that no Fermi-LAT data for PSR J2021+4026 is available between MJD 58,190 and 58,240, and we exclude this large observational gap when we bin the Fermi data. For each time bin, we describe the frequency evolution as f(t) = f(t₀) + $\dot{f}$(t₀)(t−t₀), where f, $\dot{f}$, and t₀ represent the spin frequency, the frequency time derivative and reference time, respectively. For reference time (t₀), we choose the middle point of each 70 day time bin. We perform an H-test (de Jager & Büsching 2010) to search the best frequency f(t₀) and the time derivative $\dot{f}({t}_{0})$ in the range of [3.767 Hz, 3.769 Hz] and of [−9 × 10⁻¹³ Hz s⁻¹, −7 × 10⁻¹³ Hz s⁻¹], respectively, within which the best solutions are expected from the previous analyses (Al13; Zhao et al. 2017). The bottom panel of Figure 1 shows the evolution of $\dot{f}$, and the uncertainty of each data point denotes the Fourier resolution in the periodicity search. In each data set, we can clearly identify the pulse peak in the folded light curve with the local ephemeris.

In Figure 1, we can see a change of the spindown rate ($\dot{f}$) around the epoch of the gamma-ray flux drop. No recovery of the spindown rate and gamma-ray flux are observed since ∼MJD 58,150. This indicates that PSR J2021+4026 changed the state from an HGF/LSD state to a LGF/HSD state around MJD 58,150 (2018 February 1). The averaged spindown rates at the four stages (dashed lines in the bottom panel of Figure 1), are $\dot{f}=-(7.72\pm 0.05)\times {10}^{-13}\,\mathrm{Hz}\,{{\rm{s}}}^{-1}$, −(8.15 ± 0.04) × 10⁻¹³ Hz s⁻¹, −(7.77 ± 0.05) × 10⁻¹³ Hz s⁻¹, and −(8.01 ± 0.07) × 10⁻¹³ Hz s⁻¹, which are calculated from the data points at MJD 54,800–55,820 (HGF/LSD), 55,880–56,825 (LGF/HSD), 57,100–58,120 (HGF/LSD), and 58,270–58,540 (LGF/HSD), respectively. We find that the spindown rate in the new HSD state after the 2018 state change is smaller than that of 2011 state change; the increase of the spindown rate at 2018 state change, $\bigtriangleup \dot{f}$/$\dot{f}$ = 3.1 ± 1.1%, is also smaller than the 2011 event (5.6 ± 0.9%).

Figure 2 is a close-up of the evolution of the gamma-ray flux with 10 day time bin size (top panels) and the spindown rate (bottom panels) around two state change events. The spindown 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 share the data of 60 days. For the previous event, we can see an abrupt change of the gamma-ray flux and the spindown rate near MJD 55,850 with a timescale less than 10 days (Al13). For the 2018 event, the state change would have occurred at around MJD 58,150 with a timescale similar to the previous event, although abrupt change of the new event is less clear compared to the previous event.

Zoom In Zoom Out Reset image size

The state change of PSR J2021+4026 could be triggered by the glitch, which is an abrupt change in the spin frequency of the neutron star. Figure 3 summarizes the evolution of the spin frequency around the 2018 state change. The dashed line and solid line in the upper panel of the figure show the liner functions with the slope, $\dot{f}$, in the middle column (before state change) and in the last column (after state change), respectively, in Table 1. The bottom panel shows the difference between the observed spin frequency and the predicted one from the ephemeris before the 2018 state change (dashed line in upper panel). We can see a change of the spindown rate at around MJD 58,150. With the current uncertainty given by the data, however, we cannot confirm an apparent frequency jump around MJD 58,150. If there is a glitch event, its size, $\bigtriangleup f/f$, is much less than 10⁻⁷. Our result is consistent with the 2011 event (Al13).

Zoom In Zoom Out Reset image size

Table 1.  Ephemerides of PSR J2021+4026 Derived from LAT Data of MJD 55,857–56,943, 57,062–57,565 (Zhao et al. 2017), and 58,244–58,722 (This Work)

Parameter
R.A., α	20:21:30.733
Decl., δ	+40:26:46.04
Valid MJD range	55,857–56,943	57,062–57,565	58,244–58,722
Pulse frequency, f (s−1)	3.7689669240(2)	3.7689112482(6)	3.7688306913(3)
First derivative, $\dot{f}$ (10−13 s−2)	−8.1978(1)	−7.738(1)	−8.0224(4)
Second derivative, $\ddot{f}$ (10−22 s−3)	−0.19(2)	3.5(3)	−3.1(1)
Third derivative, $\dddot{f}$ (10−29 s−4)	0.572(7)	−1.7(9)	1.4(2)
Fourth derivative, 10−37f(4) (s−4)	0.3(1)	−6(9)	⋯
Epoch zero (MJD)	56,400	57,200	58,400
rms timing residual (μs)	2199.970	2795.366	1808.901
Time system	TDB (DE405)

Note. The numbers in parentheses denote 1σ errors in the last digit.

Download table as:  ASCIITypeset image

3.2. Pulse Profile and Phase-averaged Spectrum

The previous studies (Al13; Zhao et al. 2017) fit the gamma-ray pulse profile with two or three Gaussian components, and find that the pulse shape changed at the 2011 event. In particular, the intensity ratio of Peak 1 and Peak 2 was measured with 0.54 ± 0.06 for the HGF/LSD state and with 0.24 ± 0.03 for the LGF/HSD state. The third Gaussian component was necessary in the fitting for the pulse profile of the HGF/LSD state before the 2011 event and it was not necessary for the LGF/HSD state. To examine the pulse profile after the state change in 2018, we generate the ephemeris between MJD 58,244 (2018 May 6) and 58,722 (2019 August 27; see Zhao et al. 2017 for the detailed method). Table 1 summarizes the ephemerides at different epochs; the second and third columns present the ephemerides for MJD 55,875–56,943 (previous LGF/HSD state) and 57,062–57,565 (HGF/LSD state) reported in Zhao et al. (2017), and the fourth column presents the ephemeris obtained in the current study. As shown in Section 3.1, the spindown rate after 2018 state change is smaller than that measured after the 2011 event.

Figure 4 shows a pulse profile generated with >0.1 GeV energy bands after 2018 state change and its fit to the two Gaussian model (curved line) is overlaid. The fitting parameters in the different epochs are summarized in Table 2. We find that the peak ratio in the new state after the 2018 event is P1/P2 ∼ 0.26 ± 0.12, which is consistent with the previous observation in the LGF/HSD state after 2011 state change. In the 2011 event, Peak 1 in the pulse profile changed the full width at half maximum (FWHM). In the 2018 event, we also observe the change of FWHM of Peak 1. The magnitude of the width change is consistent with the 2011 event within the error, as Table 2 shows.

Zoom In Zoom Out Reset image size

Table 2.  Parameters of Gaussian Fitting for Pulse Profiles of E > 0.1 GeV

	2011		2018
	HGF/LSDa	LGF/HSDb	HGF/LSDc	LGF/HSD
Peak 1d	0.19 ± 0.02	0.13 ± 0.02	0.195 ± 0.024	0.119 ± 0.027
Peak 2d	0.176 ± 0.007	0.174 ± 0.006	0.152 ± 0.010	0.161 ± 0.013
Peak 1/Peak 2e	0.54 ± 0.06	0.24 ± 0.03	0.494 ± 0.009	0.260 ± 0.118

Notes.

^aAl13, three Gaussian. ^bAl13, two Gaussian. ^cZhao et al. (2017), two Gaussian. ^dFWHM. ^eRatio of amplitude.

Download table as:  ASCIITypeset image

Figure 5 shows the phase-averaged spectra before (MJD 57,800–58,100) and after (MJD 58,250–58,550) the 2018 state change, and Table 3 summarizes the fitting parameters of the different epochs using a power law plus exponential cutoff function described by Equation (1). The obvious change in the spectral properties are observed for the 2018 event as well as the 2011 event: the cutoff energy decreased and the photon index increased after the state change. The similar change of the spectral properties was also observed in the 2011 event, and the observed flux level and the cutoff energy between the same states are consistent with each other. These results therefore suggest that gamma-ray emission of PSR J2021+4026 repeats the state change between two different states. If PSR J2021+4026 repeats the state change with the current observed time interval, the next recovery from the current LGF/HSD state to an HGF/LSD state may occur in 2021 and the next abrupt state change from the HGF/LSD state to a LGF/HSD state may happen in 2025.

Zoom In Zoom Out Reset image size

Table 3.  Parameters of Spectral Fitting

	2011		2018
Parameter	HGF/LSD	LGF/HSD	HGF/LSD	LGF/HSD
Flux (10−6 cts cm−2 s−1)	1.29 ± 0.01	1.15 ± 0.01	1.29 ± 0.01	1.12 ± 0.01
Cutoff energy (MeV)	2755 ± 74	2477 ± 77	2675 ± 135	2447 ± 110
Photon index	1.64 ± 0.011	1.66 ± 0.014	1.68 ± 0.02	1.75 ± 0.02

Note. The parameter in the 2011 event are taken from Zhao et al. (2017).

Download table as:  ASCIITypeset image

The timing and spectral analysis suggest that PSR J2021+4026 entered a new LGF/HSD state at around 2018 February 1. Although the spindown rate of the current LGF/HSD state is slightly smaller than that after the 2011 state change, the gamma-ray emission properties (spectrum and pulse profile) for the two LGF/HSD states are consistent with each other. Moreover, the properties of the emission properties and spindown rate between the previous two HGF/LSD states (i) before 2011 October and (ii) between 2014 December and 2018 February are also consistent with each other. The long-term observation of Fermi-LAT suggests that (i) PSR J2021+4026 may repeat an abrupt state change from an HGF/LSD state to an LGF/HSD state with a time interval of τ_int ∼ 6–7 yr, (ii) the abrupt state change occurs over a timescale of less than 10 days, and (iii) PSR J2021+4026 may stay at an LGF/HSD state for about three years and recover to an HGF/LSD over a timescale of a few months. Future continuum monitoring is required to obtain the solid conclusion of the time interval and timescale of the state change of PSR J2021+4026.

In this study, we cannot confirm the existence of the glitch at the state change. As we will discuss in this section, however, we expect that the state change is related to the neutron star's quake, which is probably related to the glitch of some pulsars. More than 180 glitching pulsars have been confirmed (Espinoza et al. 2011)¹¹ and about 50 Fermi-LAT pulsars show the glitch event (Kerr et al. 2015). However, the state changes of the gamma-ray emissions from other pulsars have not been confirmed yet. In the radio bands, on the other hand, it has been suggested that various observed phenomena, including emission mode change, nulling, and some timing noise, are related to the state change of the global magnetosphere. Lyne et al. (2010) report timing noise for 17 radio emitting pulsars, including the glitching pulsars, and they suggest that those timing behaviors are a result of switching (suddenly for some pulsars) between two different spindown rates. They also demonstrate that the change of the spindown rate of those "switching pulsars" is quasiperiodic with a timescale of days to years. Moreover, the timing noise of six pulsars is correlated to the change of the radio pulse profile, unlike the usual glitch in which no significant change of the pulse profile is observed. These results suggest that the structure of global magnetosphere of the state switching pulsars changes periodically with time. Although confirmation of the next state change from the current LGF/HSD state to HGF/LSD state is necessary, PSR J2021+4026 may also periodically switch between two different spindown states.

It has been known that the change of the spindown rate $| \bigtriangleup \dot{f}|$ of the glitching pulsars and switching pulsars (Lyne et al. 2010; Espinoza et al. 2011) is tightly correlated with the spindown rate $| \dot{f}|$ (Figure 6). This would indicate that the change of the spindown rate of switching pulsars is caused by a process similar to the glitch of some pulsars. We can see in Figure 6 that $| \bigtriangleup \dot{f}|$ of the switching pulsars tends to have larger value than the averaged values of the glitching pulsars, and the magnitude of $| \bigtriangleup \dot{f}|$ of PSR J2021+4026 is also larger than the values of the glitching pulsars with a similar spindown rate. PSR J2021+4026 therefore may provide another opportunity to study the mechanisms of the state switching.

Zoom In Zoom Out Reset image size

Besides the usual glitch process, several mechanisms have been proposed to change the spindown rate of the pulsars. For example, Ng et al. (2016) attribute the permanent-like increase in the spindown rate of PSR J2021+4026 to the change of the inclination angle by $\bigtriangleup \alpha \sim 5^\circ$. The change of the spindown rate can also be caused by the change of the global current structure. The spindown power may be described as L_sd ∼ IV_a, where I is the current and ${V}_{a}\sim 2{\pi }^{2}\mu /{P}_{s}^{2}{c}^{2}$, with μ being the magnetic dipole moment of the neutron star, is the available potential drop of the pulsar. The observed large change in the spindown rate of PSR J2021+4026 may imply the change in the global current rather than the potential drop and the current in HSD sate increases by $\bigtriangleup I/I\sim \bigtriangleup \dot{f}/\dot{f}\sim 0.03$. In the pulsar magnetosphere, the magnitude of the current would be controlled by microphysical processes (e.g., pair-creation) of the particle acceleration regions, where the electric field along the magnetic field line exists. Several possible regions have been suggested for the acceleration/emission region, e.g., polar cap accelerator (Ruderman & Sutherland 1975; Daugherty & Harding 1996), slot gap (Harding et al. 2008; Harding & Kalapotharakos 2015), outer gap accelerator (Cheng et al. 1986; Takata et al. 2016), and current sheet region (Spitkovsky 2006). Since those different acceleration regions will be connected by the global current flow (Shibata 1995), a state change of one acceleration region will affect other acceleration/emission regions and hence the global magnetosphere structure. Since it has been thought that the gamma-ray emission occurs around the light cylinder, we may expect that some physical process around the light cylinder (e.g., magnetic reconnection at the Y-point, magnetosphere instability) triggers the state change of the pulsar. However, Fermi-LAT has already discovered >200 gamma-ray emitting pulsars (The Fermi-LAT collaboration 2019), but only PSR J2021+4026 shows the switching behavior in both spindown rate and emission properties (two Fermi-LAT pulsars, PSR J2043+2740 and J0742-2822, are probably switching pulsars, Lyne et al. 2010). Since PSR J2021+4026 has typical spindown parameters in the gamma-ray pulsars, we should expect more variable gamma-ray pulsars, if the physical process around the light cylinder could trigger the state switching of PSR J2021+4026.

The state switching of PSR J2021+4026 is likely caused by change of the magnetic field structure of the polar cap region. Since the the polar cap region covers only about 0.05 (P_s/0.1 s)⁻¹% of the stellar surface, the probability of the disturbance of the polar cap region by a glitch process (e.g., crust cracking) may be low. This could explain why switching pulsar, especially switching gamma-ray pulsar is rare. In Ruderman (1991a, 1991b), the crust magnetic flux tube penetrates the type II super-conducting core. The neutron vortex lines in the core move as the pulsar spins down, and pull the magnetic flux tubes in the core with them. If the magnetic flux at the crust is immobile, the shear stress will be built up at the crust. The crust cracking will occur if the shear stress exceeds the crust elastic stress (Cheng & Dai 1998),

$$\begin{eqnarray}&&{B}_{c}\delta B/8\pi \gt {\mu }_{s}{\theta }_{\max }\displaystyle \frac{{\ell }}{{R}_{{NS}}}\end{eqnarray} \tag{ 2 }$$

where B_c is the average magnetic field strength in the quantized flux tubes at the core, δB is the magnetic field strength creating the shear stress, μ_s is the shear modulus, θ_max is the strain before the cracking, ℓ is the size of the crust, and R_NS is the radius of the neutron star. Hence, the magnetic field strength δB just before the cracking is estimated to be

$$\begin{eqnarray}&&\begin{array}{l}\delta B\sim 3\times {10}^{11}\ {\rm{G}}\ \left(\displaystyle \frac{{\mu }_{s}}{{10}^{29}\ {\rm{g}}\,{\mathrm{cm}}^{-3}}\right)\ \left(\displaystyle \frac{{\theta }_{\max }}{{10}^{-3}}\right)\\ \ \ \qquad \times \ {\left(\displaystyle \frac{{B}_{c}}{{10}^{15}{\rm{G}}}\right)}^{-1}\left(\displaystyle \frac{{\ell }}{{10}^{5}\ \mathrm{cm}}\right).\end{array}\end{eqnarray} \tag{ 3 }$$

After the cracking, the cracked plate will move to reduce the stress. The displacement (δℓ) and the time interval (τ_int) between two successive abrupt state changes may be related with the stellar radius and characteristic age (τ_sd) of the pulsar by τ_int/δℓ ∼ τ_sd/R_NS (Ruderman 1991a, 1991b). The displacement is rewritten as

$$\begin{eqnarray}&&\delta {\ell }={10}^{2}\,\mathrm{cm}\left(\displaystyle \frac{{\tau }_{\mathrm{int}}}{7\ \mathrm{yr}}\right){\left(\displaystyle \frac{{\tau }_{\mathrm{sd}}}{70\mathrm{kyr}}\right)}^{-1}.\end{eqnarray} \tag{ 4 }$$

With the parameter of PSR J2021+4026, τ_int ∼ 7 yr and τ_s ∼ 70 kyr, the displacement is estimated to be δℓ ∼ 100 cm.

With the parameters of PSR J2021+4026, the radius of the polar cap is R_pc ∼ 2.8 × 10⁴ cm, and the size of the polar cap accelerator along the magnetic field line will be of the order of ∼10³ cm (Ruderman & Sutherland 1975). Hence the sudden displacement ℓ ∼ 10² cm of the cracked plate could affect the polar cap structure. With a strong magnetic field at the polar cap region, the quantum electro dynamical processes such as the magnetic pair-creation process and/or the photon splitting process can be realized. The mean free path of such a magnetic process is sensitive to the photon energy, the strength of the magnetic field, and the angle between the propagation direction of the photon and direction of the magnetic field (Ruderman & Sutherland 1975; Timokhin & Harding 2019). If the crust cracking (or other process) would change the magnetic field structure or field strength in the polar cap region, the process affects the created current and hence the spindown rate. It is expected that the change of the current is very limited, that is, $\bigtriangleup I/I\ll 1$, since the magnitude of current I will be controlled between the particle acceleration and the screening of the accelerating field by the pair-creation process and it will be always close to the Goldreich–Julian rate (Goldreich & Julian 1969).

The quasiperiodic change of the timing residual of the radio pulsars has been widely argued within the framework of the precession of the neutron star (Stairs et al. 2000; Link & Epstein 2001; Jones 2012; Kerr et al. 2016). In the precession interpretation for a biaxial star, the precession period, P_p, is related to the spin period, P_s, by

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{s}}{{P}_{p}}=\displaystyle \frac{{I}_{0}}{{I}_{p}}\epsilon ,\end{eqnarray} \tag{ 5 }$$

where I₀ and I_p denote the characteristic moment of inertia, which neglects the small difference between different axes, of the star and of precessing body, respectively. The ellipticity is defined by $\epsilon =\bigtriangleup {I}_{p}$/I₀ with $\bigtriangleup {I}_{p}$ being the difference between the largest and smallest moments of inertia about the principal axes. With P_s = 0.265 s and P_p ∼ 6–7 yr that is observed time interval from an HGF/LSD state to an LGF/HSD state, the required ellipticity is of the order of  ∼ 10⁻⁹ if I_p ∼ I₀, and it is similar to the values for other switching pulsars. The wobble angle θ is estimated to be of the order of $\theta \sim | \,\bigtriangleup \dot{f}| /| \dot{f}| \sim 2^\circ$, which is also similar to the required angle of well studied precession candidate, PSR B1828-11 (Link & Epstein 2001).

We however argued that the change from an LSD state to an HSD state of PSR J2021+4026 happened over a timescale of ∼10 days or less, which corresponds to the phase width δΦ ∼ 10 days/6.5 yr ∼ 0.004. Such an abrupt change could be difficult to explain with a simple precession interpretation, since the variations in the spindown rate given by the precession would be expressed with the fundamental and first harmonics of the Fourier components. Jones (2012) argues that an avalanche-like process such as the pair-creation process discussed above is probably enhanced at some precession phase and magnifies the observed spindown rate. Another possibility is that a glitch of PSR J2021+4026 excites new free precession and about three years of the LGF/HSD represents damping timescale of the precession. The damping timescale is estimated as τ_d = 400−10⁴ precession period (Alpar & Sauls 1988; Jones & Andersson 2001). With τ_d = 3 yr of PSR J2021+4026, the precession period is estimated to be P_p ∼ 0.1–2.6 days. In this scenario, only a small part of the star will be precessed, I_p ≪ I₀, if the PSR J2021+4026 has a typical ellipticity  ∼ 10⁻⁸–10⁻⁹. Unfortunately, it will be difficult to confirm the precession with P_p ∼ 1 days for PSR J2021+4026 because of the radio-quiet source.

The evolution of the spindown rate is also discussed with the twisted pulsar magnetosphere, which are proposed for the evolutions of the X-ray luminosity and the spindown rate of the magnetar after the X-ray outburst (Thompson et al. 2000; Beloborodov 2009). In this scenario, the neutron star quake would twist the closed magnetic field lines and the current flowing along the twisted magnetic field lines creates a toroidal current component. The induced toroidal current can provide an additional magnetic dipole moment, μ. Since the spindown rate is proportional to μ², the twisted magnetosphere can increase the spindown rate. This scenario can explain the evolution of the X-ray luminosity after the X-ray outburst of magnetar XTE J1810-197 (Beloborodov 2009) and the high-magnetic field radio pulsar PSR J1119-6127 (Wang et al. in preparation). For PSR J1119-6127, the pulsar underwent a large glitch at 2016 July X-ray outburst (Archibald et al. 2016), and the spindown rate rapidly increased by a factor of ∼5−10 in a timescale of a month. Then the pulsar recovered to the preburst rate over the following three months (Dai et al. 2018). In the twisted pulsar magnetosphere scenario, the timescale that the twisted magnetic field exists is characterized by the evolution timescale of "j-bundle," and for a canonical pulsar it may be estimated as

$$\begin{eqnarray}&&\begin{array}{l}{t}_{t}\sim \displaystyle \frac{{\phi }_{0}{\sin }^{2}\theta \mu }{{{cR}}_{{NS}}{V}_{e}}\sim 7\ \mathrm{hr}\times {\phi }_{0}\left(\displaystyle \frac{{\sin }^{2}\theta }{0.5}\right)\\ \qquad \times \ \left(\displaystyle \frac{\mu }{5\cdot {10}^{30}\ {\rm{G}}\,{\mathrm{cm}}^{3}}\right)\ {\left(\displaystyle \frac{{{eV}}_{e}}{{10}^{12}\mathrm{eV}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 6 }$$

where ϕ₀ ∼ B_t/B is the twist angle with B_t being the toroidal field, θ is the polar angle of the j-bundle, R_NS is the neutron star radius, and V_e is the voltage to maintain the current flowing along the twisted magnetic field lines. For a canonical pulsar, the potential drop along the magnetic field line will develop so that the electrons/positrons emit the GeV photons via the curvature radiation process. The GeV photons will be converted into the electron/positron pairs by the magnetic pair-creation process and the pairs will stop the development of the potential drop along the magnetic field line. Near the stellar surface, the typical magnitude of the potential drop is of the order of eV_e ∼ 10¹² eV. With the typical dipole moment, μ ∼ 5 × 10³⁰ G cm³, of canonical pulsars, the timescale for the existence of the twisted magnetic field will be less than one day, which is difficult to explain the observed timescale ∼3 years of high spindown state of PSR J2021+4026. It is expected that if the abrupt state change of PSR J2021+4026 is caused by the twisting of the magnetosphere, the emission from the heated footprint contributes to the observed X-ray emission in an HSD state. In the previous studies of the X-ray emission from PSR J2021+4026 (Lin et al. 2013; Wang et al. 2018), there is no apparent difference between the X-ray emission properties in the HSD and LSD states. It is therefore that the twisted magnetosphere will not be the reason of the state change of PSR J2021+4026.

PSR J2021+4026 changed the state from an HGF/LSD state to an LGF/HSD state around MJD 58,150 (2018 February 1), and the state change would occur over a timescale of less than 10 days. The properties of the state change are similar to what was observed in 2011 October, although the measured spindown rate in the new state is slightly smaller than that measured after the 2011 event. The long-term evolution of the spindown rate and gamma-ray emission properties suggest that PSR J2021+4026 repeats the state change every several years. The future monitoring is required to obtain the solid conclusion of time interval of the state change of PSR J2021+4026.

More than a dozen switching pulsars have been discovered in the radio bands. PSR J2021+4026 will be the first switching pulsar discovered in the gamma-ray bands and provides a unique opportunity to understand the mechanisms of the state switching. We argued that the shear stress building up at the crust will crack the plate that covers the polar cap region. The displacement of the cracked plate, δℓ ∼ 10² cm, will affect to the structure of the polar cap accelerator, whose size along the magnetic field line will be ∼10³ cm. It would be also possible that the glitch of the neutron star excites new precession with a period of 0.1–2.6 days and the timescale of LGF/HSD state (∼3 yr) would present the damping timescale of the precession. If PSR J2021+4026 repeats the state change with the current observed time interval, the next recovery from the current LGF/HSD state to another HGF/LSD state may occur in 2021 and the next abrupt state change for an HGF/LSD state to an LGF/HSD state may happen in 2025.

J.T. and W.H.H. are supported by NSFC grants of the Chinese Government under 11573010, 11661161010, U1631103, and U1838102. A.K.H.K. is supported by the Ministry of Science and Technology of the Republic of China (Taiwan) through grants 105-2119-M-007-028-MY3 and 106-2628-M-007-005. K.S.C. is supported by GRF grant under 17302315. C.Y.H. and K.L.L. are supported by the National Research Foundation of Korea grant 2016R1A5A1013277. P.-H.T.T. is supported by NSFC through grants 11633007 and 11661161010 and U1731136.