The ORT and the uGMRT pulsar monitoring program: Pulsar timing irregularities & the Gaussian process realisation

Abstract

The spin-down law of pulsars is generally perturbed by two types of timing irregularities: glitches and timing noise. Glitches are sudden changes in the rotational frequency of pulsars, while timing noise is a discernible stochastic wandering in the phase, period, or spin-down rate of a pulsar. We present the timing results of a sample of glitching pulsars observed using the Ooty Radio Telescope (ORT) and the upgraded Giant Metrewave Radio Telescope (uGMRT). Our findings include timing noise analysis for 17 pulsars, with seven being reported for the first time. We detected five glitches in four pulsars and a glitch-like event in PSR J1825–0935. The frequency evolution of glitches in pulsars, J0742–2822 and J1740–3015, is presented for the first time. Additionally, we report timing noise results for three glitching pulsars. The timing noise was analysed separately in the pre-glitch and post-glitch regions. We observed an increase in the red noise parameters in the post-glitch regions, where exponential recovery was considered in the noise analysis. Timing noise can introduce ambiguities in the correct evaluation of glitch observations. Hence, it is important to consider timing noise in glitch analysis. We propose an innovative glitch verification approach designed to discern between a glitch and strong timing noise. The novel glitch analysis technique is also demonstrated using the observed data.

1. Introduction

Pulsars are neutron stars with exceptionally stable rotation and good timekeepers due to their huge moment of inertia, magnetic fields, and compact sizes. In particular, millisecond pulsars rival the best atomic clocks. After accounting for the steady slowdown of the underlying neutron star due to the magneto-dipole radiation or particle outflows in the form of wind, the accuracy of their spin period is even better than one part in $10^{11}$ . This allows for predicting the lengthening of the period of a given pulsar once its timing solution is obtained for a reference epoch by the pulsar timing technique. Regular spin evolution of a pulsar is occasionally interrupted by variations and rare discontinuities that cause phase deviations, the so-called timing irregularities. The long-term timing observations have revealed that the rotations of many pulsars are affected by two types of irregularities (D’Alessandro 1996): glitches and timing noise.

A glitch is a sporadic event where a pulsar shows an abrupt change in its rotation rate. Glitches in the rotation rate is usually accompanied by even larger step increases in the spin-down rate. Both of these changes tend to relax back towards the original pre-glitch values on timescales ranging from several days to a few years. In contrast, timing noise refers to a random variation in the rotational parameters of the neutron star. It is usually quasi-periodic in nature, and its frequency spectrum can be characterised by a power-law.

The first detected rotational spin-up in the Vela pulsar occurred on Modified Julian Date (MJD) 40280, February 28, 1969 (Radhakrishnan & Manchester 1969; Reichley & Downs 1969), followed by a similar event detected in the Crab pulsar in the same year on MJD 40491.80, November 1969 (Downs 1981). These events were called ‘pulsar glitches’ and motivated a deeper exploration of neutron stars’ internal structure and dynamics and shedding light on dense matter under extreme conditions. Soon after the detection of the glitches, several models were developed to explain the phenomenon; refer to reviews Haskell & Melatos (2015), Zhou et al. (2022), Antonopoulou et al. (2022) and Antonelli et al. (2022) for more details. The most successful model explaining glitches and their recoveries is based on the superfluid vortex pinning and creep mechanisms (Anderson & Itoh 1975; Alpar et al. 1984; Gügercinoğlu & Alpar 2020), which states that internal superfluid state leads to the development of a differential rotation between the observed crust and stellar interior by immobilisation of vortex lines via pinning to lattice nuclei and flux tubes in the outer core. Hence, a large amount of angular momentum gets stored. The crust spins up whenever this angular momentum is abruptly released to produce a glitch. The much slower post-glitch recovery takes place as a result of the weak interaction between the superfluid components and normal matter crust. Pulsar glitches thus probe the superfluid dynamics inside the star.

It has been more than 55 yr since the first detection, and only about 700 glitches in around 230 pulsars (combined data from the Jodrell Bank Observatory ^a Basu et al. 2022 and the ATNF pulsar catalogue ^b Manchester et al. 2005) have been recorded, making them rare events. Several types of spin variations have been observed, resulting in the classification of glitches into four categories, namely, conventional glitches, slow glitches, glitches with delayed spin-ups, and anti-glitches, for example, see the review (Zhou et al. 2022). A conventional glitch is a glitch that shows an abrupt spin-up and exhibits a typical exponential recovery followed by a linear decay of the increase in the spin-down rate. In some events, gradual increments in the frequency have been observed. These gradual change events are referred to as slow glitches (Shabanova 1998). The typical timescale for the formation of a slow glitch can vary from several months to years (Cordes & Helfand 1980). Only 31 slow glitches in seven pulsars have been reported so far (Zhou et al. 2019). Some of the glitches seen in the Crab pulsar, including its largest 2019 glitch, have displayed a unique behaviour; glitch magnitude attained its maximum value with a delay of about a few days from the beginning of the spin-up-event, forming a category called delayed-spin up glitches (Lyne et al. 2015; Gügercinoğlu & Alpar 2019; Basu et al. 2020). A sudden negative step jump in the spin frequency is called an anti-glitch (Archibald et al. 2013). Such kind of events have been identified several times, predominantly occurring in magnetars (e.g. Içdem, Baykal, & Inam 2012; Archibald et al. 2013; Şaşmaz Muş, Aydın, & Göğüş 2014; Pintore et al. 2016; Ray et al. 2019; Younes et al. 2020).

The migration and repinning of the vortices define the post-glitch behaviour. It generally exhibits an exponential recovery followed by a linear decay of the increase in the spin-down rate (Yu et al. 2013). Post-glitch recovery behaviour can be utilised to study and evaluate various statistical parameters that probe the interior structure and superfluid/superconducting dynamics inside the neutron stars. The post-glitch behaviour provides insights into the following scientific aspects: (i) the fractional Moment of Inertia of various layers that participated in a glitch, (ii) the re-coupling timescale of the crustal superfluid, (iii) the theoretical prediction of the time to the next glitch, (iv) constraining internal magnetic field configuration and equation of state of neutron stars, (v) finding a correlation between glitch observables, (vi) the range of proton effective mass in the core, and (vii) a robust internal temperature estimate for neutron stars and many more (Gügercinoğlu & Alpar 2020; Gügercinoğlu et al. 2022).

Various statistical studies have utilised the dataset of pulsar glitches to understand glitches and their correlations with various pulsar parameters. It has been observed that glitch size is at least bimodal (Eya, Urama, & Chukwude 2019; Celora et al. 2020; Arumugam & Desai 2023), categorising glitches into two groups: small and large glitches. It is not well known whether the triggering mechanisms for large and small glitches are the same. Many argue that a similar mechanism is responsible for the occurrence of small glitches and timing noise. A proportional behaviour has been observed between glitch activity and spin-down rate (Urama & Okeke 1999; Lyne, Shemar, & Smith 2000; Fuentes et al. 2017). Glitches are observed to be more dominant in younger pulsars. This is partly due to the fact that in young pulsars for which the spin-down rate is comparatively higher, the critical lag for collective vortex unpinning avalanches is easier to achieve (Gügercinoğlu & Alpar 2016).

Table 1. List of pulsars observed using ORT and uGMRT. The columns from left to right represent the pulsar Jname, pulsar period (P), the time derivative of the period ( $\dot{P}$ ), characteristic or spin-down age ( $\tau$ ), the surface magnetic field of the pulsar as inferred from the dipolar spin-down law (B), Dispersion Measure (DM), observatory name where the observation took place (Telescope), and the data availability period (in MJD).

Unlike pulsar glitches, timing noise is less understood, partly as characterising timing noise not only requires long timing baselines but can often be corrupted by noise introduced by the ionised interstellar medium (IISM), particularly in young distant pulsars, which are often found in complex environments. Nevertheless, there have been several studies to characterise the timing noise in pulsars. One such study showed that a power law adequately describes the timing noise in the millisecond and canonical pulsars (Shannon & Cordes 2010). However, for magnetars, the same power law does not hold due to the potential influence of magnetic field decay and frequently occurring crustquakes on the generation of timing noise, giving a hint that the strength of timing noise can be expressed as a combination of the rotational frequency of the pulsar and its derivatives. It was observed that the strength of timing noise differs significantly among canonical pulsars, millisecond pulsars, and magnetars, with variations spanning over eight orders of magnitude (Parthasarathy et al. 2019). It is conjectured that there is a correlation between the strength of timing noise and the rotational frequency of the pulsar and its derivatives (Shannon & Cordes 2010; Parthasarathy et al. 2020; Lower et al. 2020). Thus, detailed studies on these aspects are needed to understand the relationship between timing noise with the various properties of the pulsar, like rotation, magnetic field, age, etc.

It has been a mystery whether all types of reported glitches are real glitches or a manifestation of timing noise. The traditional glitch analysis does not consider the possible irregularities that can occur because of timing noise, which can masquerade as a small glitch. The advancements in timing noise modelling in recent years, primarily driven by Gravitational wave analysis, motivated us to develop an innovative glitch analysis methodology that is capable of differentiating a glitch from a manifestation of the timing noise. Previously, a semi-automated approach for glitch detection with a hidden Markov model (HMM) has been developed (Melatos et al. 2020) as an alternative to the traditional glitch analysis. This method has been utilised for glitch detection and defining the upper limit on the glitch size (Dunn et al. 2022b; Dunn et al. 2023a). The HMM incorporates both glitches and timing noise; however, HMM models timing noise as white noise in the torque derivative (an approximation). A majority of pulsars with glitches show timing noise with power-law spectra, which motivates including a stationery timing noise model derived from long-term timing baselines to evaluate if a small glitch is real. This is particularly important when there are large data gaps in data, making phase tracking a challenging task, or there is no apparent recovery, distinguishing a glitch from timing noise.

In this paper, we present the results of our monitoring programme using the Upgraded Giant Metrewave Radio Telescope (uGMRT) (Gupta et al. 2017) and the Ooty Radio Telescope (ORT) (Swarup et al. 1971). We detect several glitches in our data, obtain timing noise models for pulsars in our monitoring programme, and evaluate these glitches against the timing noise model using a different approach developed by us. The outline of the paper is as follows: We describe the pulsar sample and observations using the uGMRT and the ORT in Section 2. In Section 3, we explain our analysis methods. Section 4 introduces the novel glitch verification methodology. The glitch detection report, results of timing noise analysis, and novel glitch analysis demonstration are discussed in Section 5. The conclusions and prospects of future work are given in Section 6.

2. Pulsar sample and observations

We regularly monitor a sample of 24 pulsars using the ORT and the uGMRT to investigate timing irregularities in pulsars. Our sample of pulsars is listed in Table 1, consisting of pulsars of different ages (listed in column 4 of Table 1) that have displayed several glitches of various amplitudes in the past. We started with a smaller sample of frequent glitching pulsars, following the selection criteria described in Basu et al. (2020). Later, we increased our sample to include pulsars that exhibit glitches of varying amplitudes. This selection criterion enables us to test the occurrence of small glitches and to distinguish timing noise from real glitches.

We use India’s two largest radio telescopes to observe the listed pulsars. The first telescope used was the ORT (Swarup et al. 1971), a cylindrical paraboloid with an aperture of 530 m $\times$ 30 m, that is, 530 m long in the north-south direction and 30 m wide in the east-west direction, covering a declination range of –60 $^{\circ}$ to +60 $^{\circ}$ . It was built on a hill with the same latitude as the geographic latitude (11 $^{\circ}$ ), making it an equatorially mounted telescope and allowing it to track sources for about 8 hours. The telescope has a total collecting area of around 16 000 m $^2$ and an effective collecting area of around 8 600 m $^2$ . It provides an angular resolution of 2.3 deg (in right ascension) and 5.5 s (declination) arcminute for the total power system. The primary reflector is made up of a network of around 1 500 stainless steel wires running along the north-south direction. The cylindrical reflector directs the radio signals to 1 056 dipoles located at the focal line in the north-south direction. Each set of 48 dipoles is arranged into a subarray called a module. The telescope has 22 modules, which are phased to form a module beam. ORT has a beam-forming network consisting of 12 beams. Observations in the ORT are conducted at a central frequency of 326.5 MHz. The Pulsar Ooty Radio Telescope New Digital Efficient Receiver, PONDER (Naidu et al. 2015), which supports four modes for observations, provides us with real-time coherent dedispersed time series data. Dedispersion is performed with respect to the highest frequency in the band, which is 334.5 MHz. The DSPSR software (van Straten & Bailes 2011) was used for folding the raw data using the ephemeris created in our initial observations. It generates PSRFITS (Hotan, van Straten, & Manchester 2004) files. The data generally have 128 phase bins and at least three sub-integrations, used to provide reliable preliminary confirmations of glitches, even with one post-glitch epoch. Also, using three or more sub-integrations is beneficial for automated glitch detections (Singha et al. 2021a). At the ORT, bright glitching pulsars with low Dispersion Measure (DM), preferably less than 100 pc cm⁻³, are being observed with a cadence of 1–14 days.

The second telescope used is the uGMRT (Swarup et al. 1991; Ananthakrishnan 1995), an array of thirty 45-m fully steerable parabolic antennas extended over a baseline of 25 km, covering a declination range of -53 $^{\circ}$ to +90 $^{\circ}$ . The uGMRT (Gupta et al. 2017) supports observations in four bands: (i) Band 2: 125–250 MHz, (ii) Band 3: 250–500 MHz, (iii) Band 4: 550–850 MHz, and (vi) Band 5: 1 000–1 460 MHz. We have utilised Band 3, Band 4, and Band 5 for our pulsar observations. The phased array mode with an antenna configuration, including all the central square antennas and the first arm antennas, was used for observations. The data was recorded with 1 024 and 2 048 channels over a bandwidth of 200 and 400 MHz, respectively. To reduce the raw data, we have used the PINTA pipeline (Susobhanan et al. 2021), which produces PSRFITS files. The sources with high DM experience scatter broadening in the ORT and hence are chosen to be observed with the uGMRT. Additional criteria are that the glitching pulsars of variable glitch amplitudes and optimum observation time, in our case, less than 20 min each, are selected for observations. The cadence of the observations at uGMRT is 15–30 days. A few sources are common in both uGMRT and ORT either because they show exceptional post-glitch behaviour or are known to show radiative changes correlating with glitches, which is hard to show with the ORT, as it is a single polarisation telescope and can be better shown with uGMRT. We prefer to observe pulsars with low DM at lower frequencies, where they require less integration time, as pulsars are intrinsically bright at lower radio frequencies. On the other hand, the pulsars with high DM were observed using higher frequencies, primarily at the uGMRT, to mitigate the impact of scatter broadening at lower frequencies. The data typically have 128 phase bins and at least three sub-integrations. At least three sub-integrations were used to provide reliable preliminary confirmations of glitches, even with one post-glitch epoch.

3. Analysis

3.1 Timing analysis

The data obtained from observations using the uGMRT and the ORT were in PSRFITS format. Initial analysis was performed using PSRCHIVE (Hotan et al. 2004; van Straten, Demorest, & Oslowski 2012). The psrstat and psrplot functions were used for information extraction and visualisation purposes. The RFI mitigation, including removal of bad subints, channels, or bins, was executed using pazi. The pam command was used to collapse the archives in frequency and/or in time. The command paas was utilised to generate a noise-free template profile from high signal-to-noise ratio profiles observed during observations. The template profile was cross-correlated with all other observed profiles using the pat command, which allows a frequency domain cross-correlation method to obtain Time of Arrivals (ToAs). A timing model, which contains the pulsar’s intrinsic spin evolution (spin frequency and its time derivatives), astrometric terms (position and proper motion), parallax, DM, and frame-of-reference terms, is utilised in the generation of the timing solution. A good timing solution yields white and minimised residuals obtained by fitting ToAs according to the timing model with a weighted least-squares algorithm (Joshi et al. 2018). The pulsar timing software TEMPO2 (Hobbs, Edwards, & Manchester 2006; Edwards, Hobbs, & Manchester 2006) was used to perform the timing analysis of the pulsars. This traditional timing technique is highly successful, can be referred to as the foundation of observational studies in pulsar astronomy, and is widely used for various topics, including the detection of glitches. However, it does not consider the effects of timing noise in glitch verification methodology. In the traditional glitch analysis technique, glitch identification is conducted through visual inspection of the residuals, and a glitch event is defined by an abrupt positive change in spin frequency, $\Delta\nu$ , or a negative change in frequency spin-down rate, $\Delta\dot{\nu}$ . These two rapid changes simultaneously distinguish glitches from timing noise (Espinoza et al. 2014). In case of small glitches, many times, no sudden discrete negative change in $\Delta\dot{\nu}$ is observed. Hence, it is necessary to take timing noise into account to verify if such glitches are real or a manifestation of timing noise while analysing the residuals. We present a technique to estimate the timing noise and to differentiate glitches from timing noise in Sections 3.3 and 4, respectively.

The rotational evolution of a pulsar during a glitch can be expressed as a Taylor series expansion (Gügercinoğlu et al. 2022):

(1) \begin{equation}\begin{aligned}[c] \nu(t) = & \nu_0 + \dot{\nu}_0(t-t_0) + \frac{1}{2} \ddot{\nu}_0(t-t_0)^2 + \Delta \nu_g \\ & + \Delta\dot{\nu}_g(t-t_g) + \Delta \nu_d \exp(\!-(t-t_g)/\tau_d) \, \end{aligned}\end{equation}

where $\nu_0=1/P$ is the spin frequency, and $\dot{\nu}_0$ and $\ddot{\nu}_0$ are the first and second-time derivatives of frequency at epoch $t_0$ , respectively. The change in the spin frequency and its first-time derivative due to a glitch at epoch $t_g$ are represented by $\Delta \nu_g$ and $\Delta \dot{\nu}_g$ , respectively. $\Delta \nu_d$ and $\tau_d$ are the amplitude and decay timescale of the exponential relaxation part of a glitch. The procedure followed to obtain the parameters of the model, given in Equation (1), is described in Section 3.5. The long-term spin evolution for one of the pulsars in our sample, including several glitches, is given in Fig. 1.

Figure 1. The rotational evolution of the Vela pulsar (PSR J0835–4510). The top panel represents the spin evolution. The middle panel shows the frequency residuals $\Delta \nu$ , estimated by subtracting the pre-glitch timing solution, and the bottom panel represents the evolution of the spin-down rate. The vertical lines indicate the glitch epoch.

3.2 Bayesian analysis

Bayesian analysis is a method of statistical analysis that uses Bayes’ theorem to update and revise the probability for the hypothesis of an event based on new evidence or data. Bayes’ theorem is the basis of Bayesian statistics, and the foundation of Bayes’ theorem is conditional probability. We provide a succinct description of the methods used, and a more exhaustive discussion of Bayesian regression and model comparison techniques can be found in Trotta (2008), Sharma (2017), Kerscher & Weller (2019) and Krishak & Desai (2020) and references therein. The posterior probability distribution of a set of parameters $\theta$ for data D and a hypothesis/model H is given as:

(2) \begin{equation} P(\theta|D,H) = \frac{P(D|\theta,H)P(\theta|H)}{P(D|H)},\end{equation}

where

• • $P(\theta|H)$ is the prior, that is, pre-knowledge or the probability distribution of parameter $\theta$ assuming a model.

• • $P(D|H)$ is the Bayesian evidence – the probability of data or a normalisation constant,

• • $P(D|\theta,H)$ is called the likelihood – how likely the data is given the model or parameters and,

• • $P(\theta|D,H)$ is the Posterior, that is, the probability that the hypothesis is true for a given data.

We have utilised the power of Bayesian analysis for model comparison and parameter estimation of various timing noise factors, which are described in detail in Section 3.3. To select the most probable model, one has to estimate Bayesian Evidence, which is defined as:

(3) \begin{equation} P(D|H) = \int P(D|\theta,H)P(\theta|H) d\theta. \ \end{equation}

To decide the probable hypothesis between two models (say $H_{A}$ and $H_{B}$ ), we calculate the Bayes Factor (BF), which is defined as the ratio of evidence and can be written as:

(4) \begin{equation} BF_{AB} = \frac{P(D|H_{A})}{P(D|H_{B})} = \frac{\int P(D|\theta_{A},H_{A})P(\theta_{A}|H_{A}) d\theta_{A} \ }{\int P(D|\theta_{B},H_{B})P(\theta_{B}|H_{B}) d\theta_{B} \ }.\end{equation}

Jeffrey’s scale (Trotta 2008) is used to examine the strength of evidence in favour of one hypothesis or model over another. The scale classifies the evidence into various levels, ranging from ‘negative evidence’ to ‘decisive evidence’. The interpretation of the scale is as follows: An evidence ratio less than 1 indicates negative evidence. A BF between 1 and 10 infers no substantial evidence. A BF between 10 and 100 suggests strong to very strong evidence, and if the BF is greater than 100, it indicates decisive evidence in favour of one hypothesis over another. We have utilised Jeffrey’s scale for the model selection. Complex models are given preference only if they have decisive evidence. However, if the evidence is not decisive, then we followed Occam’s razor and selected the simplest model with fewer free parameters. We utilised Bayesian analysis to determine the timing noise parameters, discussed in the next section.

3.3 Timing noise analysis

The irregularities in the spin create a time-correlated stochastic function in the TOAs, known as achromatic red noise, red spin noise, or timing noise. The timing noise in our sample of pulsars has been modelled as a Gaussian stationary random process using Fourier basis functions (Lentati et al. 2013; van Haasteren & Vallisneri 2014). Unlike white noise, which has a constant power spectral density across all frequencies, red noise exhibits a power spectrum where the power decreases with increasing frequency. This results in a correlation structure where fluctuations at nearby time points are more strongly correlated than those further apart. The power spectral density of timing noise (van Haasteren & Levin 2013; van Haasteren & Vallisneri 2014; Lentati et al. 2016; Parthasarathy et al. 2019; Lower et al. 2020) for a given observing frequency f can be written as:

(5) \begin{equation} \mathcal{P}(f) = \frac{A_\mathrm{red}^2}{12\pi^2} \left(\frac{f}{f_{\text{yr}}}\right)^{-\gamma} \text{yr}^3\,, \end{equation}

where $A_\mathrm{red}$ is the red-noise amplitude, in units of $\mathrm{yr}^{3/2}$ , $\gamma$ is the power law index, and $f_\mathrm{yr}$ = 1.0/yr.

Our timing noise analysis model differs from Melatos et al. (2020) in terms that they have incorporated timing noise as a red noise in the spin-down rate. Other sources of uncertainties, such as radiometer noise and pulse jitter, can result in white noise features in the residuals (van Haasteren & Vallisneri 2014). These can be accounted by scaling the ToA uncertainties as given below:

(6) \begin{equation} \sigma^2 = F^2 (\sigma_r^2 + \sigma_Q^2) . \end{equation}

Here, $\sigma_r$ is the ToA uncertainty provided by ‘pat’ in PSRCHIVE and F is known as EFAC. The EFAC represents the radiometer noise, which is affected by Radio Frequency Interference (RFI), etc, and hence, its value can vary. However, EFAC is preferably close to unity. The term $\sigma_Q$ , known as EQUAD, is the error added in quadrature, constant for all the data points, and describes the stochastic fluctuation in the pulse profile, which can contribute additional noise to the timing data (Osłowski et al., 2011; van Haasteren & Vallisneri 2014; Kikunaga et al. 2024). While EFAC could be epoch-dependent due to RFI or variations in telescope parameters, EQUAD represents a time-independent white noise required for accurate error estimation. In timing noise analysis, we use four different combinations of noise models as listed below:

1. 1. F: Model consisting of only the white noise scaling factor (EFAC).

2. 2. WN: Model consisting of the white noise scaling, EFAC, and the white noise error in quadrature (EQUAD).

3. 3. F+RN: Model consisting of the white noise scaling factor (EFAC) along with the intrinsic pulsar red noise.

4. 4. WN+RN: Model consisting of the white noise scaling factor (EFAC), white noise error in quadrature (EQUAD), and the intrinsic pulsar red noise.

Our noise analysis is based on the Bayesian analysis method described in Section 3.2. The prior used for our analysis are presented in Table 2. The prior for red noise parameters are chosen based on those used previously (Parthasarathy et al. 2019; Lower et al. 2020) and further constrained during the analysis. The value of EFAC is observatory dependent, while EQUAD depends on the pulsar and observation frequency. The initial prior range of EFAC for uGMRT sources was chosen from previous uGMRT noise analysis (Srivastava et al. 2023) as 0.1–5. During the initial analysis, we noticed that a range from 0.1 to 3 was also suitable as a prior range for our pulsars, so we used this range for the final analysis. The same range (0.1–3) for EFAC was assumed for the sources at the ORT and was found to be a convenient choice. The initial prior range for EQUAD was chosen as -3 to -10 (Parthasarathy et al. 2019), which includes the prior ranges used before (Srivastava et al. 2023) for the uGMRT. During the initial analysis, we observed that the prior range used in Srivastava et al. (2023) and Parthasarathy et al. (2019) alone was insufficient to model the posterior distribution of several pulsars. Therefore, we expanded our range and found that the range of -1 to -10 (Lower et al. 2020) was a suitable choice for all of our uGMRT and ORT pulsars. The timing noise is modelled using Gaussian likelihood (van Haasteren et al. 2009), which is determined with Enterprise (Ellis et al. 2019) (Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE). Enterprise is a pulsar timing analysis code that can perform noise analysis, gravitational-wave searches, and timing model analysis.

Table 2. Prior ranges used for the noise analysis of our pulsars. The first column represents the parameter name, followed by the prior range and the prior distribution.

3.4 Fourier modes selection

It is important to select the optimum number of Fourier modes (M), as the noise models are sensitive to the number of Fourier bins for each pulsar (Chalumeau et al. 2022). Young pulsars experience strong red noise because of highly dynamic internal and magnetospheric activities. Therefore, it is vital to have higher modes. If the data span is $T_{span}$ (in years), the lowest frequency mode is given as 1/ $T_{span}$ , and the highest mode is M/ $T_{span}$ . Therefore, in Fourier mode selection, the sample space is 1/ $T_{span}$ , 2/ $T_{span}$ , $\ldots$ , M/ $T_{span}$ . The value of modes can be decided by following the Nyquist criteria. Ideally, for daily cadence, the highest frequency is given as $2/$ day. Therefore, the value of the highest number of Fourier modes for uniform daily observation is 730. However, our data are not uniform and has gaps between the observations. Additionally, the cadence is not strictly 1 day; for the uGMRT, it is 15–30 days, and for the ORT, 1–14 days. Accordingly, the choice for the highest mode frequency 48/ $T_{span}$ signifies a once-a-weak frequency, while the lower component with a frequency of 3/ $T_{span}$ , four months is an adequate choice for Fourier modes sample space. In principle, we should use all the modes (from 1 to M) and select the model with the largest Bayes factor. However, it can be noted that the changes in the Bayes factor for very nearby modes may not be significant. Furthermore, it may be computationally expensive to use all the modes. Similar to the previous study on modes selection (Srivastava et al. 2023), we are using a set of five Fourier modes defined as $M = 3, 6, 12, 24, 48$ (1/ $T_{span}$ ), where $T_{span}$ is the time span of the data in years for the optimum mode selection. The Fourier modes selection has been done simultaneously with the model selection. We estimated the Bayesian evidence for all the models with our chosen five sets of Fourier modes. The number of Fourier modes has been selected for GP realisation through a visual inspection of observed data and the realisation plot. The GP realisation includes observations near the glitch epochs; hence, a large number of modes, say around 100, is preferable for obtaining accurate realisation plots. However, having a very high value of components does not add any extra advantage.

Table 3. The parameters of all glitches presented in this work. The J name of the pulsars, the epoch of the glitch, the pre-glitch rotation frequency, and the rotation frequency derivative at the glitch epoch are listed in the first four columns, respectively. The last two columns present the fractional change in the rotational frequency and its time derivative, respectively. Pulsar with $*$ in JName represents the glitch-like event.

3.5 Parameter estimation

The glitch parameters have been obtained by the traditional timing analysis using TEMPO2. First, a pre-glitch timing solution was obtained, considering only the pre-glitch ToAs. Likewise, a post-glitch timing solution was obtained. Then, the glitch epoch is determined by identifying the x-coordinate of the intersection point (in MJD) of the fitted straight line in the pre-glitch and straight/quadratic curve in post-glitch residuals. Finally, the pre-glitch and post-glitch solutions were obtained with the glitch epoch as the reference. This provides the ratio of change in the rotation frequency corresponding to the post-glitch and pre-glitch solution to the pre-glitch rotational frequency. This fractional change in rotational frequency characterises as the glitch amplitude. We utilised Enterprise to model the noise in pulsars, and DYNESTY (A Dynamic Nested Sampling Package for Estimating Bayesian Posteriors and Evidences) (Speagle 2020) has been utilised to estimate log evidence for the model and Fourier mode selection. After obtaining the preferred model and an optimum number of modes, we used PTMCMCSampler (Ellis & van Haasteren 2017) for sampling and parameter estimation. The Python library Corner (Foreman-Mackey 2016) has been utilised to plot the posteriors.

4. The Gaussian process realisation

Are all glitches real, or are some of them manifestations of timing noise? How to differentiate between a real glitch and a glitch that appeared because of strong timing noise? We have developed a novel approach to answer these questions. It is a five-step process:

1. 1. Use the traditional glitch analysis process, that is, using the pulsar timing technique to detect glitch-like events.

2. 2. Model the timing noise as a red noise process separately for the pre-glitch and post-glitch sections using the Bayesian inference.

3. 3. Utilise the results of Bayesian noise analysis to reconstruct a GP realisation for pre-glitch and post-glitch parts.

4. 4. Subtract the median GP realisation signals from the observed signals, which should result in residuals distributed as white noise.

5. 5. Use the whitened residuals to verify if the glitch is real or not. A real glitch is characterised by a jump soon after the glitch epoch in the whitened residuals.

To plot the GP realisation curve, each value of red noise amplitude $A_\mathrm{red}$ and spectral index $\gamma$ , that is, chains obtained from the timing noise analysis, are used. The Gaussian process is the family of curves obtained in the timing noise analysis, and the median value of the posterior for $A_{red}$ and $\gamma$ will give the median curve; if one subtracts the median curve from the residuals, then we obtain whitened residuals. For our analysis, we utilised the Python LA-FORGE (Shafiq Hazboun 2020) package for the GP realisation.

The demonstration of the process is given in Section 5.3, where we used a large glitch event (in PSR J0835–4510) and a small glitch (in pulsar J1847–0402), which appears like a glitch, but it turned out to be a manifestation of the timing noise. The GP realisation helped us to verify these events. A step increase in the rotational frequency due to strong timing noise or lack of phase connection (due to large observation gaps) can lead to a small glitch. Hence, we recommend utilising this technique to test all glitches.

5. Results and discussions

We started the glitch monitoring programme using the ORT and the uGMRT in 2014. The first results of this monitoring programme, with 11 glitches in 8 pulsars, were reported in Basu et al. (2020). Here, we present the detection of five glitches in four pulsars. Additionally, we present timing noise results for 20 pulsars, including three glitching pulsars, where the timing noise has been estimated for the pre-glitch and the post-glitch regions. Furthermore, we demonstrate our novel method to differentiate glitches from strong timing noise using the GP realisation.

5.1 Glitches

We report the detection of 5 glitches in our 4 pulsars and a glitch-like event in PSR J1825–0935. The results are tabulated in Table 3. The time evolution of two glitches, in PSR J0742–2822 and PSR J1740–3015, is being reported for the first time. Some of the glitches have also been reported by us earlier as Astronomer’s Telegram ^c (Singha et al. 2021b; Grover et al. 2022; Grover et al. 2023) as well as by other pulsar-monitoring programmes (Basu et al. 2022; Lower et al. 2021; Zubieta et al. 2024). Most of the glitches reported here are large glitches. In addition to these large glitches, we detected around 8 glitch-like events. However, our new analysis technique implied that seven of these are unlikely to be glitches as these appear to be consistent with timing noise in these pulsars. In Section 5.3, we critically examine a small glitch in PSR J1847–0402 with our new technique and show that this glitch is consistent with timing noise. We find that often such an event looks like a glitch either due to low cadence observation or big gaps between the observations. Therefore, we decided not to report these as glitches. We are only reporting one glitch-like event in PSR J1825–0935, which appears to be a slow glitch. More description of the event is provided in Section 5.1.5. Now, we concisely describe glitches in each pulsar.

5.1.1 PSR J0534+2200

The Crab pulsar (PSR J0534+2200) is one of the frequently glitching pulsars, with about 30 glitches observed since its first glitch. This pulsar is monitored regularly with a cadence of 1–3 days using the ORT and a biweekly cadence using the uGMRT. In our monitoring programme, we have detected multiple glitches, including the largest ever glitch seen in this pulsar (Basu et al. 2020). In this work, we present the glitch parameters of the 30th glitch reported for the Crab pulsar (Shaw et al. 2021), estimated to occur at MJD 58686.4(8). We calculated the fractional increase in the spin to be $26.6(3)\times10^{-9}$ , and the time derivative of spin to be $0.6(1)\times10^{-3}$ . The timing residuals and time evolution of the frequency and its derivative are presented in Fig. 2. The previous glitches from our monitoring programme, which were published in the previous paper (Basu et al. 2020) are not included in the present work.

Figure 2. A glitch seen in PSR J0534+2200 on MJD 58686. The top panel shows the residuals. The middle panel displays the $\Delta \nu$ evolution, and the bottom panel presents the evolution of $\Delta \dot\nu$ . The glitch epoch is shown by the vertical line.

5.1.2 PSR J0742–2822

With an age of 157 Kyr, PSR J0742–2822 is an intriguing old glitching pulsar, which has displayed many glitches of various amplitudes. It is also known to show frequent changes in the magnetospheric state that may occur due to a glitch (Keith, Shannon, & Johnston 2013), making it an excellent candidate to investigate. It is currently monitored at the ORT with a cadence of 1–3 days. We observed the largest glitch detected in this pulsar on MJD 59839.8(5), which was detected by the Automated Glitch Detection Pipeline implemented at the ORT (Singha et al. 2021a) and reported by us in the Astronomer’s Telegram (Grover et al. 2022). The glitch was also reported by Shaw et al. (2022a), further confirmed by Dunn et al. (2022a), and Zubieta et al. (2022b) on the Astronomer’s Telegram. This article presents a more complete analysis of all our data and final updates to our preliminary results reported in the Astronomer’s Telegram (Grover et al. 2022). The fractional rise in the spin is observed as $4\,299(10)\times10^{-9}$ , and the time derivative of spin is estimated to be $90(8)\times10^{-3}$ . We also note around 35-day recovery in the glitch, which is being reported for the first time and matches with recent reports (Zubieta et al. 2024). The glitch parameters are broadly consistent with preliminary estimates reported by other groups in the Astronomer’s Telegram(Shaw et al. 2022a; Dunn et al. 2022a; Zubieta et al. 2022b), while the timing residuals, frequency, and frequency derivative evolution, shown in Fig. 3, are being reported for the first time.

Figure 3. A glitch seen in PSR J0742–2822. The top panel shows the residuals. The middle and the bottom panel show the evolution of $\Delta \nu$ and $\Delta \dot\nu$ . The vertical line indicates the glitch epoch.

5.1.3 PSR J0835–4510

The Vela Pulsar, PSR J0835–4510, is a fantastic pulsar for glitch observations, as it shows regular glitches of similar amplitudes (large glitches with fractional sizes $\gtrsim10^{-6}$ ). We have been regularly monitoring this pulsar using the ORT with a cadence of 1–3 days and a biweekly cadence using the uGMRT. We present two recent glitches observed in the Vela pulsar using the ORT. Both these glitches are large glitches, which are the characteristic features of the glitches seen in this pulsar. The timing residuals, the evolution of the frequency, and its time derivatives for both the glitches are shown in Fig. 4. We detected its twenty-third (G23) (Kerr 2019; Lower et al. 2020; Gügercinoğlu et al. 2022) and twenty-fourth (G24) glitches (Sosa-Fiscella et al. 2021; Zubieta et al. 2023). We estimated the glitch epoch for G23 as MJD 58517(7), the fractional change in frequency is $2\,471(6)\times10^{-9}$ , and the fractional increase in the time derivative of frequency is $6(2)\times10^{-3}$ ; however, the reported values in the literature for the fractional increases are slightly higher. This is probably due to the number of pre-glitch and post-glitch epochs used for estimations and uncertainty due to slight differences in epoch estimation. For glitch G24, the glitch epoch is estimated as MJD 58417.6(1), the fractional change in frequency is $1\,235(5)\times10^{-9}$ , and the time derivative of frequency is $8.0(7)\times10^{-3}$ . The fractional change in frequency is again slightly lower than the value reported on the Jodrell Bank Observatory catalogue; again, this could be because of the number of pre-glitch and post-glitch epochs used for estimations and uncertainty due to slight differences in epoch estimation.

Figure 4. Two glitches observed in PSR J0835–4510 on MJD 58517 and 59417, respectively. The top panels represent the residuals. The middle and the bottom panels display the evolution of $\Delta \nu$ and $\Delta \dot\nu$ . The vertical line indicates the glitch epoch.

Figure 5. A glitch detected in PSR J1740–3015. The top panel shows the residuals. The middle panel displays the $\Delta \nu$ evolution, and the bottom panel presents the evolution of $\Delta \dot\nu$ . The vertical line represents the glitch epoch

Figure 6. A glitch seen in PSR J1825–0935. The top panel presents the residuals. The middle panel displays the $\Delta \nu$ evolution. The vertical line indicates the glitch epoch.

5.1.4 PSR J1740–3015

PSR J1740–3015 is one of the most frequently glitching pulsars with 37 glitches. It is presently in the second position in glitch count, which makes it a remarkable candidate for studying timing irregularities. It is currently monitored at uGMRT with a cadence of 15–20 days. We report the detection of a large glitch in this pulsar observed on MJD 59 934.8(5), which was initially announced in the Astronomer’s Telegram by Zubieta et al. (2022a) and subsequently confirmed in the Astronomer’s Telegram by Dunn et al. (2023b) and Grover et al. (2023). We update our preliminary results reported on Astronomer’s Telegram (Grover et al. 2023) by augmenting subsequent observations. We calculated the fractional increase in the spin as $327(1)\times10^{-9}$ , and the time derivative of spin is $1.3(3)\times10^{-3}$ . While we believe these are more accurate estimated parameters, they are broadly in line with the parameters reported earlier. The timing residuals, spin, and spin derivative evolution are given in Fig. 5 and are being presented for the first time for this new glitch.

5.1.5 Glitch-like event in J1825–0935

PSR J1825–0935 is one of the few pulsars which shows slow glitches (Zou et al. 2004; Shabanova 1998; Shabanova 2005; Zhou et al. 2019) along with conventional glitches. Slow glitches are characterised by an extended rise in their frequency evolution, with a typical rise timescale ranging from months to years. We have been monitoring this pulsar regularly using the uGMRT and the ORT, and we detected a slow glitch-like event also reported in Lower et al. (2020) and Singha et al. (2022). The residual and spin evolution plots are given in Fig. 6.

Figure 7. Timing residuals plots for 17 pulsars used for single pulsar noise analysis presented in this work.

To date, only $\sim$ 31 slow glitches have been reported (Zhou et al. 2019) and the mechanisms behind slow glitches are not yet fully understood. However, there are several theoretical models that suggest the occurrence of slow glitches. One such model suggests that the cause of slow glitches is the oscillation between two phases of an anisotropic superfluid (Peng, Liu, & Chou 2022). Others proposed that a slow glitch occurs due to vortex bending oscillation in the spin-down rate, generated either by magnetospheric changes or crustquakes (Gügercinoğlu, Köksal, & Güver 2023). Slow glitches are dominantly observed in old pulsars with low spin-down rates. As in such systems, superfluid recoupling timescale would be longer, and this, in turn, affects the migration rate of vortex lines once they become unpinned. When transported radially inward driven by a crustquake, completion of the angular momentum transfer process and achieving the equilibrium configuration for vortex lines take a longer time due to the old age of a pulsar (Gügercinoğlu & Alpar 2019). The occurrence of slow glitches in PSR B0919+06 concurrent with cyclic changes in its spin-down rate is a firm demonstration (Shabanova 2010).

PSR J1825–0935 is a very interesting pulsar because it also exhibits profile changes concurrent with its glitches. It is known to display interpulse, drifting subpulses, and mode switching between two emission states: Burst/B-mode and Quiet/Q-mode (Fowler, Wright, & Morris 1981; Morris, Graham, & Bartel 1981; Fowler & Wright 1982; Gil et al. 1994). In the burst mode, an additional faint precursor unit adjacent to the intense main pulse is visible, whereas in the quiet mode, the precursor unit exhibits minimal emission, and the interpulse appears brighter. The timing study of this pulsar has complications because of mode switching. According to Lyne et al. (2010), the slow glitches are actually not related to the glitch phenomenon and are consequences of the switching between the two magnetospheric spin-down modes. The observed oscillations in the spin-down rate, correlated with pulse profile changes or induced after glitches (Shaw et al. 2022b), support the view that either magnetospheric processes or crustquakes lead to a change in the strength of coupling of the interior superfluid torque acting on the neutron star surface (Gügercinoğlu & Alpar 2017; Gügercinoğlu et al. 2023).

It has also been suggested that slow glitches are not unique and are a manifestation of the timing noise (Hobbs, Lyne, & Kramer 2010). However, if slow glitches were manifestations of merely timing noise in pulsars, they should not be rare (Zhou et al. 2019). Our novel glitch verification methodology is capable of differentiating a real glitch from timing noise. However, the observed data using ORT and uGRMT do not have sufficient pre-glitch observations for timing noise analysis. Hence, the GP realisation has not been utilised for this pulsar. Additionally, due to the complications in the timing from mode switching, PSR J1825–0935 is not a good example to answer if slow glitches are real. A rigorous analysis accounting for all the reported slow glitches is required to answer whether they are real glitches or a special case of timing noise.

5.2 Timing noise

We now present the noise analysis of 20 pulsars in our sample, including 3 glitching pulsars. The timing residuals of 17 pulsars, with no glitch observed using the uGMRT and the ORT are shown in Fig. 7. The pulsar timing package, Enterprise, was used to obtain the timing noise parameters, and PTMCMCSampler has been used for sampling. A description of the four different combinations of noise models and various model parameters is given in Section 3.3. We performed the Bayesian model selection using these models. In Table 4, we present the $\ln$ (BF) of the various models with respect to the simplest model in terms of the number of parameters. The best model was selected for each pulsar based on the values of the Bayes Factor. The bold value of the number against a model indicates that it is the most preferred model. This preferred model was selected based on the values of the BF and for its simplicity. In cases where the BF of two or more models are similar, we selected the model containing the least number of parameters (Occam’s razor). The preferred models and their corresponding parameters for the 20 pulsars are tabulated in Tables 5 and 6. The posterior distribution for the parameters of the most preferred noise models for our pulsars are given in 7. The posterior distributions are plotted with 68% and 95% credible intervals for all the cases.

Table 4. The values of $\ln(\text{BF})$ with respect to the simplest model, in terms of parameters (F) for the 20 pulsars in our sample. Bold values against a model indicate it is the preferred model for the corresponding pulsar. This preferred model has been selected based on the values of the BF and the number of free parameters in the model. The pulsars marked with $*$ have experienced glitches. The model comparison results correspond to the pre-glitch region, and for the post-glitch region, the preferred model remains the same; however, the value of $\ln(\text{BF})$ may vary.

Three of our pulsars, namely, PSR J0525+1115, J1720–1633 and J1919+0021 show evidence of only white noise (EFAC + EQUAD) in the timing residuals. The following pulsars, J0358+5413, J0729–1448, J0729–1836, J0922+0638, J1532+2745, J1731–4744, J1740–3015, J1847–0402, J1909+0007, J1910–0309, J2022+2854, J2219+4754, J2346–0609 show evidence for red noise and do not prefer the inclusion of EQUAD in the noise model. On the other hand, PSRs J0528+2200, J0742–2822, J0835–4510, and J1825–0935 preferred the full white noise (EFAC + EQUAD) along with the red noise model. Some pulsars prefer the inclusion of EQUAD while others do not because, for some pulsars, the observed timing residuals contain extra noise contributions from epoch to epoch pulse shape changes or jitter (Osłowski et al. 2011) that are not well-characterised by the formal measurement uncertainties. Adding an EQUAD term helps to account for this extra noise, leading to more accurate parameter estimates and uncertainties. The noise results for PSRs J0525+1115, J0729–1836, J0742–2822, J0835–4510, J0922+0638, J1720–1633, J1731–4744, J1740–3015, J1825–0935, J1847–0402, J1909+0007, J1910–0309, J1919+0021, J2346–0609 have been reported in literature (Lower et al. 2020; Parthasarathy et al. 2019).

The white noise model for PSRs J0525+1115, J1720–1633, J1919+0021 and the red noise model for PSRs J0729–1836, J0922+0638, J1731–4744, J1825–0935, J1847–0402, J1909+0007, J2346–0609 were also preferred in Lower et al. (2020). The only difference was in pulsar J1910–0309, which preferred the white noise model in their results and showed evidence of red noise in ours. This is probably due to the availability of a larger data span in our case. The reported red noise parameters for most of the pulsars, for example, J0729–1836, J0922+0638, J1731–4744, J1825–0935, J1847–0402, were within the error bars of the parameters reported by them. However, the timing noise parameters are different for a few pulsars, such as PSRs J1909+0007 and J2346–0609. This inconsistency can be due to the unmodelled interstellar medium effects in both ours and Lower et al. (2020) analyses. The timing noise (or red noise) parameters for pulsars J0358+5413, J0528+2200, J0729–1448, J1532+2745, J1910–0309, J2022+2854, and J2219+4754 were not reported. Therefore, we present the new timing noise results for at least seven pulsars.

Table 5. The noise parameters obtained for our sample of pulsars that have not experienced any glitch. The 1st column represents the pulsar J name, followed by the time span of observation, the 4th column represents the most preferred model for the pulsar, the 5th and 6th columns present the white noise parameters, and the 7th and 8th columns are the red noise parameters. The last column corresponds to the number of glitch-like events observed in the respective pulsar, which are unlikely to be real glitches as they appear consistent with timing noise in these pulsars.

Table 6. The noise parameters obtained for our sample of pulsars that have experienced glitches. The 1st column represents the pulsar J name, followed by the time span of observation, the 4th column represents the most preferred model for the pulsar, the 5th and 6th columns present the white noise parameters, and the 7th and 8th columns are the red noise parameters. The last column corresponds to the number of glitch-like events observed in the respective pulsar, which are unlikely to be real glitches as they appear consistent with timing noise in these pulsars.

We have also analysed the timing noise in glitching pulsars in the pre-glitch and the post-glitch regions for three pulsars for the first time to the best of our knowledge. The preferred models and their corresponding parameters for three glitching pulsars are listed in Table 6, and the timing noise posteriors of the most preferred noise models are given in Appendix A (Figures 11–12). A significant variation in timing noise parameters has been observed before and after the MJD 59839.8 glitch in PSR J0742–2822 and the MJD 59417.6 glitch in PSR J0835–4510. This is probably due to the consideration of the exponential post-glitch recovery region in the noise analysis. However, as no exponential recovery was present for the MJD 58517 glitch in PSR J0835–4510, the noise results in the pre-glitch and the post-glitch regions are consistent.

Figure 8. Verification of glitch using our new technique employing GP realisation for PSR J0835–4510.

Figure 9. The plot for PSR J0835–4510 residuals obtained after subtracting the GP realisation from the timing residuals revealing the clear signature of a glitch.

Figure 10. Verification of glitch using our new technique employing GP realisation for PSR J1847–0402.

5.3 Differentiating glitches from timing noise : GP realisation

As mentioned in Section 4, we devise a novel method to distinguish real glitches from timing noise by taking into account the timing noise in glitching pulsars. We illustrate this process here with two extreme examples.

5.3.1 PSR J0835–4510

We begin with the glitch in the Vela Pulsar, which is a large glitch; hence, there is a high probability that it is a real glitch. Additionally, we have a high cadence ( $\sim$ 1–3 days) data for this glitch, as it is being observed at the ORT regularly. We started with the pre-glitch timing solution for this pulsar and utilised it to investigate the timing noise and obtain timing noise parameters/chains. Further, timing noise results were used to produce a GP realisation. The results are illustrated in Fig. 8a. The top panel of the figure displays the observed data points alongside the median GP realisation. The bottom plot illustrates the subtracted data obtained by subtracting the GP realisation from the observed data points. The residual data appears consistent with zero mean white noise. Now, in a glitch event, we anticipate observing a jump in the post-glitch region after eliminating the timing noise. Fig. 8b presents the GP realisation for the post-glitch region. The curvature observed in the bottom plot is attributed to the residual data being fitted in tempo2. This process generated a negative compensated component to counteract the glitch jump event, resulting in the curvature observed in the white noise signal. The curvature gives an initial hint that the glitch will probably be real. If we exclude the initial few weeks of observations, primarily representing the exponential recovery phase, and focus on the linear recovery region, we anticipate observing a flat line once more, as shown in Fig. 8c. Finally, we extrapolate the residuals described in Fig. 8c and subtract the extrapolated value for the initial observations or the exponential recovery phase. This yields the final glitch plot, providing evidence that the glitch is real. The final glitch plot is shown in Fig. 9.

5.3.2 PSR J1847–0402

Here, we discuss the case of PSR J1847–0402. Three small glitches (glitch amplitude $\sim 10^{-10}$ ) have previously been reported for this pulsar. We detected a glitch-like event on MJD 59883(6) using timing analysis with a glitch amplitude of $0.18(3)\times 10^{-9}$ , which is close to the range reported previously. We utilised the GP realisation to verify this glitch. The results are shown in Fig. 10, authenticating that this event was not a real glitch but a manifestation of the timing noise. It is possible that all other small glitches reported in this pulsar appeared because of strong timing noise. Hence, it is important to utilise our new method of employing GP realisations to investigate small glitches, which are expected to appear frequently because of strong timing noise. Overall, we detected around 8 glitch-like events, seven of which appeared like glitches because of large observation gaps or strong timing noise, and the 8th one is a glitch-like event in J1825–0935, as discussed in Section 5.1.5.

6. Conclusions and future work

Pulsar glitches and timing noise probe into the dynamics of superfluid interiors and help us to understand the physical conditions prevailing inside neutron stars. We present the recent results of our glitch monitoring programme using the uGMRT and the ORT. This includes the detection of five glitches in 4 pulsars. The time evolution of two of these glitches has been presented for the first time, showing significant recovery from the glitch. Most of the glitches presented in this work are large glitches. Additionally, we report noise analysis for 20 pulsars using Bayesian analysis techniques given in Tables 5 and 6. The timing noise analyses for 14 of the pulsars from our sample were presented before (Lower et al. 2020). We are reporting new results for noise analysis of 6 pulsars and red noise results for J1909–0309 for the first time. Furthermore, we have investigated the timing noise in the pulsars that have shown glitches and presented the timing noise parameter values in pre-glitch and post-glitch regions. A significant variation in the pre-glitch and post-glitch red noise parameters has been observed (given in Table 6) for the cases where exponential recovery has been considered for the noise analysis.

The continuous increment in the detected glitches may shed light on many correlations that are not well understood yet, for example, correlation with inter-glitch time, and may provide more information about the mechanism responsible for the occurrence of glitches. The reported glitch data sets are utilised to study statistical properties like differentiating the glitching pulsar population from the non-glitching population. Hence, it is important to report real glitches only. A discontinuity in observations or strong timing noise may manifest itself as a small glitch and result in a false detection, which should be taken care of. To tackle this problem, we have suggested a novel glitch verification method that can differentiate between a real glitch and a pseudo-glitch, which is a manifestation of any other phenomenon.

About $\sim$ 35% of the total number of glitches reported have relatively small amplitudes. This suggests that many of the small glitches that have been reported could be due to strong timing noise in pulsars. Hence, it is very important to use our updated glitch verification methodology to make sure that any considered glitch is real and not a manifestation of timing noise.

The proposed GP realisation technique for glitch analysis is capable of differentiating between a real glitch and a pseudo-glitch, which is demonstrated in Section 5.3 with basic noise models, consist combinations of white noise and achromatic red noise. These simple noise models serve as a good starting point. However, in future endeavours, it may be necessary to employ more comprehensive timing noise models by including several effects, such as the effect of the interstellar medium, which can play a significant role for some of the pulsars, for example, the Crab pulsar.