Singularity Detection on Forbush Decrease at High Latitude Stations During Geomagnetic Disturbances

We analyzed the behavior of Cosmic Ray (CR) intensity during five geomagnetic events over 9-day periods of varying disturbance, using ground-based CR measurements from the Bartol Research Institute neutron monitor program that occurred on 25 October – 2 November 2003, 23rd – 31 st July 2004, 13 th – 21 st May 2005, 13 th – 21 st Mar 2015, 19 th – 27 th June 2015. A nine-day time frame choice gives us a reasonable time interval for data analysis before, during, and after the main event(s). These events have been deliberately chosen, intending to track Forbush Decrease (FD) at high latitude stations following geomagnetic storms of different intensity and duration from perspectives of their origin and geoeffectiveness. FD was observed when the magnetic fields entangled in and around CME exert a deflecting effect on galactic cosmic radiation, resulting in a sudden reduction in their intensities. The results revealed that the FD plunged between -4% and -17% in the chosen events. The FD examined was abnormal, and a multi-stage decrement in FD was observed during the event period. Furthermore, we have used the Discrete Wavelet Transform (DWT) technique to detect singularity on CR intensity at the stations described. The first three decomposition levels have proved sufficient to isolate the transients in CR intensity in conjunction with different nature and intensity events, ranging from intense to super intense geomagnetic storms to High‐intensity long‐duration continuous auroral electrojet activity (HILDCAA) event. In addition, it has been found that the percentage decrease in cosmic ray intensity was correlated with the Dst index during the process, as indicated by the cross-correlation technique. No noticeable lag has been found between the parameters discussed, which indicates a good association between the IMF Bz and the Dst index and the FD. Possible physical explanations supporting the results are discussed.


Introduction
Cosmic rays (CRs) are eminently energetic charged particles that mainly comprise ~89 % of protons, ~10% of helium nuclei, and ~1 % denser elements (Mewaldt, 2010;Singh et al., 2011). Cosmic rays are categorized into different types based on their origin and interaction. Some are spawned from the Sun, known as solar energetic particles (SEP), while others originate as galactic cosmic rays (GCR) from astrophysical origins outside the solar system (Mathpal et al., 2018). The modulation of cosmic ray intensity (CRI) is influenced by various solar and interplanetary variables (Usoskin et al., 2008;Belov, 2008). Huge magnetic disturbances reach the Earth in the form of either interplanetary coronal mass ejections (ICMEs) or solar flares, deviating incoming cosmic rays and decreasing their intensity (Cane et al., 1996;Heber et al., 1999;Cane, 2000). Such events act as interplanetary shock drivers. CRIs substantially drop with the emergence of these disturbances (Mathpal et al., 2018). These reductions in the CRIs were first observed by Scott. E. Forbush and termed as Forbush decrease (FD) (Forbush, 1938;Mathpal et al., 2018). This fluctuation in intensity of cosmic rays is observed due to the interaction of CRs with the solar wind and the frigid magnetic field enclosed inside it. Lockwood (1971), Cane (2000), and Richardson (2004) address FD events thoroughly in detailed reviews.
Forbush decrease is characterized as the transient decrease in cosmic ray intensity caused by the disturbances in interplanetary conditions emanating from coronal mass ejections (CMEs), solar flares, and high-speed solar wind streams (Cane, 1995;Joselyn, 1986;Singh et al., 2008). Various researchers have depicted two types of Forbush decreases as non-recurrent decreases and recurrent decreases (Lockwood, 1971;Cane, 2000). Non-recurrent decreases are associated with the ejection of solar flares and/or CME, while recurrent decreases are associated with the coronal holes related to high-speed solar wind streams (Lockwood, 1971;Singh et al., 2008). Although both FDs and geomagnetic disturbance are manifestations of a given interplanetary disturbance, their magnitudes are not always proportional to each other (Kane, 1977;Singh et al., 2008). FDs are described as a sudden onset (often with a complicated time structure) that reaches a minimum within a day, followed by a more gradual recovery period that can last anywhere from a few days to a few weeks (Cane, 2000;Zhao and Zhang, 2016). The earliest phase of FD is observed prior to the arrival of the interplanetary disturbance at Earth as a pre-increase of CR intensity caused by the acceleration of high energy charged particles on the outer boundary of an ICME and a pre-decrease due to magnetic connection between the Earth and the region of FD inside of ICME (Belov, 2008). When the propagation of an interplanetary disturbance is faster, and the magnetic field is higher, the Forbush decrease is faster (Belov, 2008). The majority of FDs are sporadic and caused by ICMEs. In this case, CRI decrease is created by the expansion of a disturbed solar wind region that is partly screened from outside by strong and/or transverse magnetic fields (Belov et al., 1997;Belov et al., 1999). Finally, when an ICME propagates beyond the Earth orbit, we see the phase of Forbush effect (FE) recovery, showing that an expanded disturbed region continues to modulate CRs (Belov, 2008). In other words, as the shock travels away from Earth, its influence on CR modulation at 1 AU weakens with heliocentric distance, allowing the cosmic ray intensity to gradually recover (Le Roux and Potgieter, 1991).
The FDs characterize the investigation of transient CR modulation, space weather forecasts, and physical processes relating to energetic particle circulation through the perturbed medium (Kuwabara et al., 2009;Zhao and Zhang, 2016). On all time scales, cosmic-rays provide a wealth of information about processes occurring in the universe in general and in the vicinity of the Earth in particular (Belov, 2000;Cane, 2000). For a long time, we have been aware of the benefit of FE observation for space weather tasks, but effective tools for practical use of FDs information are not developed so far. Without FE observation, the picture of solar and heliospheric storms would not be complete and transparent. Solar wind, magnetosphere, and cosmic ray disturbances are all connected and induced by the same active processes on the Sun. As a result, investigating them together is reasonable and beneficial. Furthermore, understanding which physical processes are most critical for particle transport includes the analysis of FDs.
In this work, we seek to identify the footprints of Forbush decreases (FDs) in the CR flux. The cause and general shape of Forbush decreases of cosmic rays are relatively well understood; however, the in-depth analysis of FDs in the frequency-domain is essential for the knowledge of frequency content in the provided signal. So, the purpose of this work is to validate the wavelet technique as a way to pinpoint the FD regions during geomagnetic disturbances. DWT is used as an alternative to the Fourier transform, as Xu et al. (2008) illustrated. Furthermore, DWT decomposes time-series into both frequency and space domains (Klausner et al., 2016). This work is presented as follows. Section 2 describes the data sets and methodology; section 3 illustrates the results and discussion, and section 4 concludes the results.

Data and Methodology
We considered five interplanetary origin geomagnetic intensity events, including intense and super-intense geomagnetic storms and a HILDCAA event. The types and the cause of space weather events considered in the study are listed in Table 1 It is well established that the IMF-Bz plays a significant role in determining the strength of the solar wind -magnetosphere coupling since it regulates the fraction of the energy in the solar wind flow extracted by the magnetosphere. When Bz is strongly negative, magnetic reconnection between the IMF and the geomagnetic field creates open field lines, allowing mass, energy, and momentum to be transferred from the solar wind to the Earth's magnetosphere (Dungey, 1961). Many aspects of the Earth's magnetosphere and high-latitude ionosphere are known to be dependent on Bz in a frame of reference that enables the orientation of the Earth's magnetic axis to be determined (e.g., GSM frame) (Davis et al., 1997). The Dst index is a well-known geomagnetic index for describing the symmetric equatorial ring current. The Dst index is calculated on a one-hour basis using data from four geomagnetic observatories located at middle to low latitudes (Du et al., 2010). Thus, IMF-Bz indicates the variation on a geomagnetically disturbed event, whereas the Dst index relates to its strength and occurrence. According to Gonzalez et al. (1994), geomagnetic storms can be classified as: weak (−50 < Dst ≤ −30 nT), moderate (−100 < Dst ≤ −50 nT), intense (−250 < Dst ≤ −100 nT), and very intense (Dst ≤ −250 nT). HILDCAAs are the geomagnetic activities associated with the periodic southward interplanetary magnetic field (IMF) components of Alfvén waves, which satisfy the conditions: (a) The AE index peaks to 1,000 nT at least once during the event, (ii) AE does not drop below 200 nT for more than 2 hr at a time, (iii) the conditions last for minimum two days, and (iv) they occur outside the main phases of geomagnetic storms (Tsurutani & Gonzalez, 1987;Khanal et al., 2019). Several studies have shown that these events have profound consequences for the magnetospheric/ionospheric environment (Hajra et al., 2015;Chapagain et al., 2012;. Therefore, it is of our interest to study the onset effect of the geomagnetic disturbances of different intensities on the CR intensity. A comparative analysis method was applied in the study, i.e., our study is based on observing FDs and validating the performance of wavelet transform in each type of event. We took hourly cosmic ray intensity data of five different high latitude stations, corrected for pressure and efficiency from the Bartol Research Institute Neutron Monitor Database (http://neutronm.bartol.udel.edu). The neutron monitor stations considered in the study are listed in Table 2, with the last column indicating the cutoff rigidity of each station. In the study, we have considered pressure-corrected cosmic ray intensity data in the form of percentage change from the average value. The global location of the neutron monitor stations is pointed in the global map shown in Fig. 1. We used this cosmic ray intensity data to apply the DWT Technique for tracking the FD regions observed by high latitude stations.

Methodologies
The Discrete Wavelet Transform (DWT) Technique is used here to analyze the periods of magnetic disturbances. This technique pinpoints the singularities in the spotted transition regions in the cosmic ray during geomagnetic storms. As high-frequency signals are procured at higher latitudes displaying the singularity patterns, this study mainly highlights the stations at the higher latitudes. The cross-correlation method has also been used to depict the comparison of the geomagnetic index, Dst, and the cosmic ray intensity during the five selected storms (Table 1).

Wavelet Transform
We often encounter two kinds of signals in the physical world: stationary signals and nonstationary signals. The former is a signal whose frequency content does not change over time, while the latter is a signal whose frequency content changes over time. Most raw signals are described in their time domain, and the graph representing them is simply a time amplitude representation of that signal. The analysis of signals is required to understand any physical mechanism, and analyzing any signal involves accessing information about the signal's frequency content and the time at which that frequency occurs. Various mathematical transformations have been developed as tools to process these signals and obtain information about them. Fourier transformation (FT) is one of them which deals with the frequency content of signals (Gao and Yan, 2011;Strang, 1993). It does not, however, provide any detail on when these frequencies appear. Thus, unless one is just interested in understanding the spectral component, Fourier transformation is not a suitable technique for studying the characteristics of non-stationary signals. Wavelet analysis has been surging as an inclusive technique to scrutinize the non-stationary signals and transform the data, operators, or functions into distinct frequency or scale components (Strang and Nguyen, 1996;Foufoula-Georgiou and Kumar, 1995;Daubechies, 1992;Ruskai et al., 1992;Chui, 1992. This approach can work with the dynamic window; a window will be narrowed automatically to observe the highfrequency content and be widened to capture the high-frequency content in the signals (Strang and Nguyen, 1996;Foufoula-Georgiou and Kumar, 1995). During the analysis, the mother wavelet is decomposed into series of basis functions containing dilated and translated mother wavelet functions. In Figure 2, we can see the time and frequency resolutions of the different transformations. During the signal decomposition, it is represented in a time-frequency plane as shown in fig 2, and this plane is layered with cells, called Heisenberg cells or boxes whose area is controlled by the Heisenberg uncertainty principle (Meyer, 1990;Daubechies, 1992;Sifuzzaman et al., 2009). This principle dictates that one cannot measure with arbitrarily high resolution in both time and frequency (Ruskai et al., 1992;Daubechies, 1992). Thus, the initial time series has a high timedomain resolution but no frequency domain resolution. This implies that we can discern very small features in the time domain but not in the frequency domain. Opposite to that is the Fourier Transform, which has a high resolution in the frequency domain and zero resolution in the time domain Domingues et al., 2005). The Short Time Fourier Transform has a medium-sized resolution in both the frequency and time domain (Chui, 1992;Chui, 1994). This shows WT is far better than the FT and STFT for time-frequency analysis but within the Heisenberg condition (Chui, 1992). This shows that as the box's width is made narrower, the box's height is elongated, which illustrates that the high-frequency component can be well localized in time, but the frequency localization is lost as frequency increases (Chui, 1992). Also, as the width of the box is made broader, the height of the box is getting narrower, and this illustrates that the low-frequency component is well localized in frequency but poorly localized in time.
The analysis of each part with a scale of analogous resolution depicts that the wavelet is very slender at high frequencies, whereas it is broad at low frequencies. Consequently, wavelet transform is a magnificent tool having time-frequency localization that contemplates momentary high frequencies phenomena. Wavelet transforms were developed by integrating the efforts of physicists, mathematicians, and engineers (Morlet, 1983;Klausner et al., 2011).

Discrete Wavelet Transform
The Discrete Wavelet Transform (DWT) is a linear multilevel efficient transform that works with data vectors with lengths that are an integral power of two. This transformation is very popular in data compression methodologies (Mallat, 1989;Daubechies, 1992;Mendes et al., 2005;Domingues et al., 2005) due to its fast algorithm and good local representation in multilevel. However, depending on the discretization scheme used, DWT may or may not have the redundant representation (Daubechies, 1992;Mendes et al., 2005;Klausner et al., 2011). To avoid redundancies, select wavelet functions that form an orthogonal basis and define the DWT, which employs discrete scale (j) and localization (k) values as input parameters (Puccio et al., 1994;Meyer, 1990): The wavelet functions in these sets are orthogonal, and their respective functions are translated and dilated.
In wavelet analysis, signals f(t) is represented by series like where = (2 − − ) are called mother wavelet (Daubechies, 1992;Chui, 1992), and the wavelet coefficients are given by, The coefficients are termed as "details," which specify the local errors and display the signal regularity of higher-frequency structures as per Klausner et al. (2014a).
The amplitude of the wavelet coefficients can be linked to abrupt signal variations or "details" of higher frequency as a property of wavelet analysis. (Meyer, 1990;Daubechies, 1992;Chui, 1992). The simplest orthogonal mother-wavelet function is the Haar wavelet. The DWT using Haar wavelet detects abrupt variations, i. e., one localization feature in the physical space . For instance, the Haar wavelet has a compact support in physical space and large support in Fourier space. In other words, Haar wavelets have a better localization on time than on frequency what is expected by Heisenberg's uncertainty principle.
In practice, the DWT is always implemented as a filter bank, a cascade of high-pass and lowpass filters (Daubechies, 1992;Mendes et al., 2005). This is because filter banks are a very efficient way of splitting a signal into several frequency sub-bands. To apply the DWT on a signal, we start with the smallest scale. As we have noticed before, small scales refer to high frequencies. This means that we start by looking at high-frequency behavior. At the second stage, the scale increases with a factor of two (the frequency decreases with a factor of two), and we analyze behavior around half of the maximum frequency. At the third stage, the scale factor is four, and we are analyzing frequency behavior around a quarter of the maximum frequency. And this goes on and on until we have reached the maximum decomposition level.
The Daubechies wavelet function of order two is chosen because it can detect the first-order local disturbances in the signal and its derivatives. As they form an orthogonal system, no redundant information is stored within the wavelet coefficient (Klausner et al., 2011). One property used here is that in the wavelet domain of smooth data, few localized structures (many wavelet coefficients with small amplitude) can be neglected, and one still has a good representation of the data . Therefore, in these cases, we have a compact representation of the data, and we can identify the disturbed regions just by looking at the amplitudes of the wavelet coefficients (Sifuzzaman et al., 2009). Thus, the discrete wavelet technique has emerged as an effective tool in the atmospheric signal analysis used to spot the minuscule particularities in the data (Klausner et al., 2011, Adhikari et al., 2017.

Cross-Correlation
Cross-correlation is a technique of studying different variables to determine similarities and draw similar relative characteristics, which can reveal new information (Liou et al., 2001;Usoro, 2015;Poudel et al., 2020;Silwal et al., 2021). Pearson's correlation coefficient (r) is the best method known to date, in which the correlation coefficient ranges from -1 to +1. The curve towards ±1 indicates a very strong correlation, whereas the curve around zero represents a moderate or less correlation (Katz, 1988).
The cross-correlation technique was used to analyze and compare the Dst index and percentage decrease in CR intensity during various geomagnetic storm events at various stations around the world. The time (hours) ranges from -50 to +50 on the horizontal plane. The time (hours) scale is dependent on the order in which the variables are used, so it can be used to determine which index leads or lags before and after they are correlated. By analyzing the cross-correlations between signals, one aims to better diagnose and understand the whole system (Yuan et al., 2015. Because of its simplicity, the traditional cross-correlation analysis (TCA) is now the most common method in statistics.

Relation of FD with solar wind parameters (IMF-Bz) and geomagnetic indices (Dst)
The dynamicity of Forbush Decrease during five major geomagnetic storms from the 23rd and 24th Solar cycle at the five above-mentioned stations is analyzed in this section. In addition, interplanetary parameters IMF-Bz and Dst -Index were studied concurrently along with Forbush decrease to describe the manifestation of geomagnetic disturbances.

Event_1: 25 October to 2 November 2003
This FD event exhibits one step in the shock-sheath region and multistep decrement during the Magnetic Cloud transition (Raghav et al., 2017). The Sun was inactive for a few days before the storm event. Active region 486 (S17W09) discharged a major flare of X10/2b at ~20:49 UT on 29 October as registered by GOES X-ray flux (1-8 A˚) data. NOAA (National Oceanic and Atmospheric Administration) also noted the intense radio emissions accompanied by the flare (Alex et al., 2006). In the upper panel of Fig. 3, daily values of cosmic ray intensity (%) deviation of the five neutron monitor stations have been plotted. In the event, cosmic ray intensity is sharply decreasing (16% -McMurdo, 15% -Newark, 17% -the South Pole, 14% -Thule, 15 % -FortSmith) at the late afternoon of 29 October in all five stations and recovery is also fast (within two days). Furthermore, it is known that during an FD, a decrease in CR intensity with amplitude inversely proportional to the cutoff rigidity related to the location of the CR station is observed (Samara et al., 2018).

Event -2: 23 June to 31 July 2004
Three geomagnetic storms occurred in July 2004 due to a sequence of coronal mass ejections suggesting the FD has one step in the shock-sheath region and multistep decrement during the Magnetic Cloud (MC) transit (Pedatella et al., 2008;Raghav et al., 2017). Three intense storms occurred in succession during 23rd -28th July. This unusual storm had the first storm, with a sudden storm commencement at 10:36 UT and the storm's main-phase onset at 19:00 UT on 22 July (Li et al., 2010).
Scrutinizing the middle panel of figure 4, Dst reached a minimum of −99 nT at 02:00 UT on 23 July. During the recovery of the first storm, the subsequent coronal mass ejection escalated the storm again. At 22:00 UT on 24 July, the second main phase was seen where Dst reached −136 nT at 16:00 UT on 25 July. The IMF Bz-component was inclined towards north from 17:00 UT on 23 July till 11:00 UT on 24 July. After that, IMF Bz swung south for an hour before sharply turning north. Over the next 9 hours, the IMF-Bz oscillated between south and north. On 24 July, at 22:00 UT, IMF-Bz turned southward and remained southward with a significant negative digression for over eighteen hours. The IMF-Bz inclined was altering between north and south for the next 4 hours and turned southward, before turning acutely northward at 03:00 UT on 26th July. The third storm started at 23:00 UT on 26 July; the main phase began at 05:00 UT on 27 July. A minimum Dst of −170 nT was observed at 13:00 UT on 27the July. At about 05:00 UT on 27 July, IMF-Bz turned southward and remained southward for about twelve hours, supporting the main phase as indicated by Dst values. This discrete phenomenon was supported by the work of Ngwira et al. (2012).
In the upper panel of figure 4, values of cosmic ray intensity (%) deviation of the five neutron monitor stations exhibit cosmic ray intensity is sharply decreasing (13% -McMurdo, 10% -Newark, 16% -the South Pole, 10% -Thule,12 % -FortSmith) during the early morning of 27 July. Drop-in cosmic ray intensity is supported by the IMF-Bz and Dst parameters on the middle and lower panel, indicating FDs. Steady recovery can be seen in all five stations in the next few days. The SSC occurred at 03:00 UT on 15 May as recognized in the Dst index (52nT) of the bottom panel of figure 5. Value of Dst reached -229 nT at 07:00 UT. At the onset of 15 May, the IMF-Bz suddenly inflated northward until 06:00 UT; as seen in the middle panel of figure 5, then it remained at the same pace. Following that, the IMF-Bz abruptly turned southward, peaking at -36.5 nT at 06:00 UT. IMF-Bz remained at this stage for about 2 hours (until 08:00 UT) before recovering. The IMF Bz reached its zero stage around 20:00 UT and continued northward, similar to the explanation (Sharma et al., 2011). The sudden decrease of IMF-Bz to −36.5 nT at 06:00 UT supports the abrupt decrease in Dst index value starting the main phase of the geomagnetic storm. (Dashora et al., 2009) found the same phenomenon during the onset of the main phase of the storm. Dst suddenly decreased to -247 nT around 08:00 UT on 15 May, indicating high geomagnetic activity.
FD was observed when the magnetic fields entangled in and around CME exert a deflecting effect on the galactic cosmic radiation causing a sudden reduction in the count rate from the neutron monitors (Chauhan et al., 2008). FD (7% -McMurdo, 7% -Newark, 9% -the South Pole, 7% -Thule, 9 % -FortSmith), as shown in the upper panel of the figure detected by neutron monitors on 15 May 2005, was an irregular event in response to intense solar flare associated with CME. A major geomagnetic storm followed the event. The onset time of this FD was ~7:00 UT, with a recovery time in the next six days, as shown in figure 5.

Event -4: 13 March to 21 March 2015
On 15 March, a magnetic filament arose jointly with a slow C9-class solar flare from sunspot AR2297 (Nava et al., 2016). Moderate 2-hour long C9 flare was released by the region assisted by a fast partial halo CME with an Earth-directed component (Tomova et al., 2017). The blast catapulted a CME into space. Earth's magnetic field was delivered a glancing blow during the late hours of 17 March as suggested by NOAA Model. The CME bashes the Earth on Saint Patrick's Day, igniting the most intense geomagnetic storm of the 24th Solar Cycle (Tomova et al., 2017).

Event -5: 19 June to 27 June 2015
The mid-June of 2015 was marked by higher activity in one of the largest active regions (AR) 12371 in the 24th solar cycle (Samara et al., 2018). The AR initiated two flares on June 21 and 22, 2015, accompanied by CMEs starting (Kravtsova et al., 2017). In addition, multiple substorms and a major geomagnetic storm occurred on Earth on 22 June 2015 (Astafyeva et al., 2018).
Two interplanetary shocks (IS) arrived at 06:00 UT and 18:00 UT due to geomagnetic disturbances on 22 June. Dst value soared to 13 nT and -8nT, respectively. The first IS caused an enhancement in substorm activity, as shown by the Dst index (-51 nT). The second IS was followed after a sudden increase in the Dst index and started to decrease rapidly. Dst value was minimum (-204 nT) at 4:00 UT on 23 June 2015, as shown in the bottom panel of figure 7. During the initial phase of the storm on the 22nd, the polarity of IMF-Bz changed several times. At 19:00 UT, the IMF-Bz minimum was -26.3 nT; after an hour, it turned north and reached +15.8 nT at 20:00 UT. After that, IMF (Bz) inclined towards northward till it reached 0.3nT at 01:00 UT on 23 June. From 01:00 UT to 05:00 UT (main phase of the storm), it was seen southward inclination and then northward from 05:00UT to 07:00 UT. Change in polarity of IMF-Bz continued till 12:00 UT of 24 June before the recovery phase started, as shown in the middle panel of figure 7.
The worldwide network of neutron monitor stations recorded a small Forbush effect with amplitude of the modulation of the CR neutron component's intensity of ~6-9% on 23 June 2015 (Kravtsova et al., 2017). Five neutron monitor stations recorded a Forbush effect (5% -McMurdo, 4% -Newark, 6% -South Pole, 5% -Thule, and 4 % -FortSmith) on 23rd June, 2015. An unusual structure of the FD on 24 June was seen, second FD was spotted less than three days after the first. After the minimal value of the CR intensity on 23 June 2015, a very short recovery phase began. It lasted from few hours to almost one day and was followed by a second FD that began at 12:00 UT on 24 June 2015. Finally, the usual recovery came after the event. Two FDs and the Dst connection can be observed in the figure. The outcome was an interesting and unusual formulation of the galactic CR flux, which appeared as FDs series (Samara et al., 2018).

DWT Analysis
The present study takes account of the Daubechies wavelet coefficient amplitude of order two of three levels. The use of this technique relies on the fact that it aids in determining changes in local regularity as well as the presence of very small and localized variations as signal discontinuities (Mallat, 1989;Mendes et al., 2005;Domingues et al., 2005;Ojeda et al., 2014). Daubechies (1990), Mallat (1997), Ojeda et al. (2014), Klausner et al. (2014b), Adhikari et al. (2017) all provide a more detailed description of the DWT analysis. We can detect first-order local disturbances in a signal and its derivatives using the wavelet technique (Adhikari and Chapagain, 2015). The disturbed periods are shown by the trough-like depression of the analyzed function, where the amplitude of the square wavelet coefficient is the local indicator of signals (Ojeda et al., 2014;Klausner et al., 2020). The signal exhibits singularity patterns of the CR flux during the period of geomagnetic disturbances, achieving maximum square wavelet coefficient (SWC) . SWC with the lowest value foreshadows the absence of FDs. Figure 8 shows the behavior of the SWC for discrete wavelet transforms related to the CR intensity on a super intense geomagnetic event that occurred on 29-31 October 2003 on five different neutron monitoring stations. Each diagram's uppermost panel depicts CR intensity, while the lower three panels, from top to bottom, depict the amplitude of squared wavelet coefficients at three decomposition levels (j=1,2,3) denoted by d1, d2, and d3, respectively. Only the local regularities detected on the three levels of decomposition could be considered due to a physical process, in this case, an SH + MC propagation, shown in Table 1, all identified by the increase of the wavelet's coefficients three decomposition levels.

Event_1: 25 October to 2 November 2003
In fig 8, it was possible to identify the discontinuities related to the FD during the geomagnetically disturbed period. One can observe that the highest relative amplitudes of the coefficients coincide during the mid-day of 29 October, showing the cosmic ray flux is globally affected, at least in the selected time resolution. The present analysis shows that for higher latitude magnetic stations, larger amplitude coefficients are more frequent during the initial and main phase of the storm in all three decomposition levels (d1, d2, and d3). The larger value of SWC accounts for the fact that negative deviation in CR intensity became more profound with the expanded disturbed region, which continues to modulate CRs. These coefficients have the ability to highlight the singularities associated with the provided signal, recorded by monitoring stations in response to geomagnetically disturbed periods .  In fig 9, it was found that the wavelet technique can localize the behavior of the FD during the intense geomagnetic storm, started at 22:49 UT on 26 July, the main phase began at 05:00 UT on 27 July. For this event, it is possible to detect the presence of FD in association with a physical process, in this case, a SIRs with ICME propagation, shown in Table 1, all identified by the increase of the wavelet's coefficients at the three decomposition levels. The highest relative amplitudes of the coefficients are coincident in time (26 July -27 July) in all stations, as witnessed from d1, d2, and d3. The abrupt variations of the cosmic ray intensity during the intense geomagnetic storm are emphasized by the larger amplitude of square wavelet coefficients.  The singularity patterns related to a rapid decrease in galactic CR intensity are identified during the HILDCAA event time intervals. As mentioned earlier, it was a rare event in response to an intense solar flare associated with halo CME, followed by a major geomagnetic storm. Following the Dst index, the SSC occurred at 02:39 UT on 15 May and the onset time of this FD was 7:00 UT. The results show that the highest relative amplitudes of the coefficients are coincident in time (15 May) in all stations. Among all these stations, the DWT analysis of Newark shows an entirely different pattern. It shows a more transient variation in the decomposition level d1 of the wavelet transforms on 16 May. This is because DWT acts as a cascade of high-pass and low-pass filters, and the first decomposition level keeps track of highfrequency details. It means that we first analyze the high-frequency behavior of the signal. Then, at the second and third decomposition levels, we analyze low-frequency behavior, which helps to capture the slowly varying changes within the signal (Daubechies, 1992, Xu et al., 2021. The higher amplitudes of the square wavelet coefficients indicate the singularities or discontinuities in records of CR intensity associated with the enduring time of the HILDCAA event Ojeda et al., 2009). Thus, the first three decomposition levels have proved to be sufficient to detect the physical discontinuities in the CR intensity during the active period of HILDCAA.  However, all the mentioned stations did not show the same pattern in the first decomposition level. According to the present analysis, the intensity of CR for FD events appears to depend on the onset time of a given FD event, which differs at widely separated NM stations. These differences are influenced by the longitudes of the NM stations rather than their latitudes (Fenton et al., 1959;Kichigin et al., 2018). However, during the most intense geomagnetic storm, all stations' overall singularity pattern was identified (Tomova et al., 2017). Furthermore, the cosmic ray flux recorded in all the stations accounted for a slight Forbush decrease (see fig 11, upper panel), and the wavelet coefficients of the second and third decomposition levels show a better time localization indeed.  In figure 12, the singularity patterns on CR intensity are identified with an intense geomagnetic storm at high latitude stations, and one can observe that the highest relative amplitudes of the coefficients are coincident at two different times, viz. (22 -23) Jun and 24 June in all stations. The maximum value of the wavelet coefficient in all three decomposition levels accounted for the first FD during 22-23 June, which was identified for the first interplanetary shock (IS).

Event_5
Similarly, the singularity patterns of CR intensity for the second FD, which occurred on 24 June was well defined by the third level of decomposition in all stations except for McMurdo. McMurdo presented the highest wavelet coefficient at the first and second decomposition levels for the second FD. This shows that the results obtained are encouraging because all three decomposition levels well localize the FD precisely at the time of occurrence of the main phase.
Furthermore, the FD events in UT and LT are classified as Simultaneous FD and Nonsimultaneous FD based on the overlap of the main phase event (Lee et al., 2015). Simultaneous FD is detected all over the world, regardless of where the NM stations are located (Oh et al., 2008;Oh & Yi, 2009). The main phase of a simultaneous FD event overlaps in UT for all widely separated NM stations, whereas the main stages occur at various times at different NM stations in LT. It implies that the difference in CRI may differ according to the local time of the NM stations.
The DWT technique has been shown to provide "zoom in" to localize the specific behavior, i.e., FDs associated with events of different nature and intensity ranging from intense geomagnetic storm to super intense geomagnetic storm HILDCAA event. We also noticed that these stations do not show a similar singularity pattern in different decomposition levels. However, we did not find a significant difference between the most relevant decomposition levels at the higher latitude stations. This behavior may be related to differences in interplanetary cause of the event, onset time of the FD event, local time, magnetic activity, and large spatial anisotropy in CR flux (Lee et al., 2013;Lockwood, 1971;Oh et al. 2008). They also proposed that a simultaneous FD event occurs when a stronger magnetic cloud passes by the Earth through the central portion of the magnetic barrier, whereas a non-simultaneous FD event occurs only when a weaker magnetic cloud passes through the magnetosphere on the dusk side. , respectively. The results of cross-correlation analysis are useful for gaining insight into the relationship between the Dst index and FDs and testing the dataset. We aligned two-time series, one of which is delayed relative to the other, as its peak occurs at the lag at which the two-time series are best correlated, that is, the lag at which they best line up (Menke and Menke, 2012).

Cross-Correlation Analysis
In Fig. 13, the curves attained a maximum cross-correlation coefficient of >0.85 at time lag 0 hours during the super intense (HILDCAA) event of 13 May -21 May 2005. Similarly, the Dst index shows a strong association with FD events with a maximum correlation coefficient of > 0.90 at relatively less time lag in the range 0 to +2 hours during the intense event of 13 -21 March 2015. Again, the cross-correlation curves of the Dst index and deviation in CR intensities for the intense event of 23 -31 July 2004 reached a maximum value in the range ~ 0.60 -0.77 at comparatively higher lags in the range +6 to + 10 than at the aforementioned events. The positive lag suggests that the second signal features (CR intensity) occur at later times than corresponding features in the first signal (Dst). The less value of lag occurs when a sudden decrease in CR intensity responds quicker to the change in Dst. On the other hand, we found correlation coefficients in the range 0.80 to 0.82 at time lag -11 to -5 hours and 0.80 to 0.90 at time lag (-3 to -1) hours, respectively during the events 25 October -2 November 2003 (super intense event) and 19 -27 June 2015 (intense event). The lag is negative when features in the second signal (CR intensity) occur earlier than corresponding features in the first signal (Dst). In all the observed cases, as visible in figure 13, Dst correlates with the deviation in cosmic ray flux at a relatively more minor lag value, concluding that all the stations exhibit a similar pattern during a particular event of study. It is interesting to notice the high correlation found in all the series, even though the datasets do not exactly match.  As the CME overwhelms the Earth, we witnessed FD due to the deflecting effects of distorted magnetic fields detected by cosmic ray detectors (Arunbabu et al., 2013). During the process, cosmic ray flux was collated with IMF-Bz values and the Dst index. However, no noticeable lag has been found between the parameters discussed for the super intense (HILDCAA) event of 13 May -21 May 2005, intense event of 13 -21 March 2015, and super intense events of 25 October -2 November 2003. It indicates a clear association between the IMF-Bz and the Dst index and the FD, leading to the view that FDs and Dst depression could be correlated with the same physical mechanism. However, the different nature of time lag during the intense event of 23 -31 July 2004 and super intense event of 25 October -2 November are probably due to the effects of many nonlinear processes and external forcing, causing the signals to have multi-scale structures and non-stationarities (Barnett, 1991;Huang et al., 1998). Overall, the neutron monitoring stations' response to the FD is reasonably rapid in these events, as reflected by a very low time lag value. This indicates that the findings obtained are encouraging because, in the near future, by incorporating precise measurements of lead or lag in the response of the Dst index with Forbush Decrease, the space weather activity related to FE could be easily predicted by placing an array of cosmic ray stations at equatorial to high latitude. Tracking of FDs is the most natural and straightforward method to investigate CME and other solar activity manifestations over a long period. Thus, FE data provide information on previous and coming storms.

Conclusion:
We inspected short-term decreases in the intensity of the Galactic cosmic rays observed from the Earth, termed as Forbus Decreases, during five major geomagnetic events from the 23 rd and 24 th Solar cycle. FD plunged in the range of -4% to -17% in the chosen events. FD reviewed were irregular and multistep decrement in FD was found during the period of the study. FD was observed when the magnetic fields entangled in and around CME exert a deflecting effect on galactic cosmic radiation, causing a sudden reduction of CR's intensity.
FDs and the Dst connection have been distinctly observed via the southward inclination of IMF-Bz during the main phase of FDs. Cosmic ray flux was collated with IMF-Bz values and Dst index during the process in all the selected events. Slight lead/lag was observed between the parameters discussed, closely matching the original time series analysis. After considering lead/lag on the particular event, correlation of IMF-Bz and Dst index with FD were studied.
Using the Daubechies orthogonal wavelet of order two of three levels, the first three decomposition levels of the wavelet components were closely examined to detect singularities in the CR intensity for the five different events ranging from intense to super intense geomagnetic storms to HILDCAA event. The preliminary remarks in this analysis can be summarized as follows: • Often, the information that cannot be readily seen in the time domain can be seen in the frequency domain. So, the wavelet technique can easily localize the behavior of FD in association with the disturbed periods of the geomagnetic storm or HILDCAA event. • We noticed a strong correlation of IMF-Bz and Dst index with FD for the five neutron monitoring stations studied. Wavelet coefficients identified it at all three levels of decomposition. The small amplitudes observed in the wavelet coefficients corresponding to cosmic ray intensity's usual behavior without geomagnetic disturbances. On the other hand, large amplitudes attest transients in cosmic ray intensity caused by the disturbances in interplanetary conditions emanating CMEs, solar flares, and high-speed solar wind streams. Thus, the first three decomposition levels have proved sufficient that the wavelet tool is useful to study singularities and fluctuations associated with CR intensity, as an indication of FDs and their influence and development during different phases of the geomagnetic disturbances. From the physical point of view, this is an indication of the real meaning of this tool.
Cross-correlation analysis presented a clear association between the Dst index and the FD. This study leads to some valuable findings, suggesting that the depression of CR intensity and Dst index may be linked with the same physical mechanism. In all the mentioned events, the neutron monitoring stations' FD response to the Dst of magnetic observatories is fairly rapid, as reflected by a meager time lag value. Overall cross-correlation and DWT results point to the fact that these approaches provide an alternative way of assessing the global effects of geomagnetic storms on CR intensity and could be used as a sophisticated tool to forecast the emergence of these unique events by tracking change in the Dst index and relative deviation in FD. However, the present study dealt with just a few geomagnetic events, five neutron monitoring stations, and one magnetic station. Further studies using more events and more stations to do a complete analysis would be needed.
The discrete wavelet transform has revealed to be an effective tool in the time localization of the cosmic ray flux using neutron monitoring station data, and it may be the most convenient to represent times series with abrupt variations or steps, i.e., very small and localized variations as the discontinuities in the provided signal used for the study.