Application of Wavelet Methods in the Investigation of Geospace Disturbances: A Review and an Evaluation of the Approach for Quantifying Wavelet Power

Abstract

Geospace disturbances refer collectively to the variations of the geomagnetic field and the trapped particle populations in the near-Earth space. These are the result of transient and recurrent solar activity, which consequently drives the variable solar wind. They may appear in multiple timescales, from sub-seconds to days, months and years. Wavelet analysis is one of the most popular, and powerful, methods in the study of these variations, as it allows for the local decomposition of non-stationary time series in frequency (or time-scale) and time simultaneously. This article is a review of the wavelet methods used in the investigation of geomagnetic field oscillations, which underlines their advantages as spectral analysis methods and demonstrates their utilization in the interdependence of multiple time-series. Lastly, the proper methodology for the accurate estimation of the power inferred from geophysical signals, applicable in quantitative studies, is included and is publicly available at the database of the University of Athens.

1. Introduction

The scales and amplitude of fluctuations in the solar-terrestrial environment have been of interest since the discovery of the solar rotation and the 11-year solar activity cycle, quantified by the sunspot (or Wolf) number or the 10.7 cm radio flux. Second (5.5 years) and third (3.7 years) harmonics of the basic 11-year periodicity have also been reported by several studies. Nevertheless, spectral analysis has revealed significant variations in solar activity on time-scales shorter than the sunspot cycle, such as:

• The Quasi-Biennial Oscillation (QBO) of ∼1.3–2 years, which is associated with double peaks of the solar cycle (see review by Hathaway [1]).
• The 154-day periodicity, which was first detected by Rieger et al. [2] in gamma-ray flare activity and has been associated with r-mode waves (Rossby waves) in the solar interior [3].
• Ephemeral periodicities due to the synodic solar rotation [4,5].

The solar wind magnetized plasma, being the expansion of the solar corona, is modulated by the fluctuations of the solar magnetic field; thus, it is affected by the quasi-periodic behaviour of the solar cycle [6,7,8,9,10]. The interaction of the solar wind with the terrestrial magnetosphere is, in turn, responsible for a similar effect of the solar cycle on the geomagnetic activity, and this chain interaction is known as solar–terrestrial coupling. Nevertheless, the geomagnetic activity is highly complex—significantly more so than expected from the solar energetic phenomena driving.

This is due to the high speed solar wind streams from coronal holes on the magnetosphere; these introduce an additional driver, out of phase with the sunspot cycle, into the Sun–Earth connection mechanism [11,12,13]. In addition to the already mentioned effects of the solar–terrestrial coupling, the periodic driving of the terrestrial magnetosphere is also apparent in cosmic ray modulation [14,15].

With the advent of wavelet analysis, which became popular in the 1990s [16], the wavelet methods have offered significant contributions in the research of the solar–terrestrial coupling and the geospace disturbances, consequential to their unique capabilities regarding data analysis; these are related to the decomposition of data into different frequency or scale components localized in time (e.g., Balasis et al. [17], Mandea and Balasis [18], Balasis and Mandea [19], Kunagu et al. [20], Zaourar et al. [21]). This article attempts an extensive review in the theory of the wavelet methods and their application in studying the quasi-periodic variations as well as the temporal behaviour for several aspects of solar–terrestrial coupling.

The manuscript is organized as follows: Section 2 begins with a historical outline of the refinement of the Fourier transform towards the wavelet time–frequency analysis including, in passing, the transitional step of the short Fourier transform (STFT) (Section 2.1). The remainder of Section 2 (Section 2.2) is a brief exposition of the methodology employed in the study of geospace disturbances, which is centred on the continuous and the cross-wavelet transform (CWT and XWT, respectively) furthered by the wavelet coherence. In Section 3, the results of the application of wavelet methods in the behaviour of the geomagnetic field, the magnetospheric particles and the Ultra-low frequency (ULF) waves are portrayed. Finally, in Section 4, a proper approach for the accurate estimation of wavelet power spectral density (PSD) is presented.

2. Theoretical Approach of Time–Frequency Analysis

2.1. Historical Context

This section is a brief exposition of the development of the wavelet area starting from eminent ancestors, such as the Fourier Transform (FT) and its shorter relative, the STFT. This presentation is somewhat qualitative as details are given below in Section 2.2. An informal but interesting historical narrative is provided by Daubechies [22] with the warning that every exposition of this type, hers not exempted, is significantly influenced by personal view and experience. Other historical accounts exist, such as Farge [23], Strang [24], Lee and Yamamoto [25] and Lau and Weng [26], which are mostly focused on the Fourier to wavelet transform transition and often appear at the introduction of research publications or tutorials. Worthy of particular mention are the more extensive presentations by Graps [27] and Akujuobi [28], Chapter 3; the latter includes a chronological list of the evolutionary stages leading to different generations of wavelets.

The transformation from the time to frequency domain is a time-honoured and widely used methodology applied, among others, in the analysis of signals, either continuous or discrete. The main advantage of this type of transforms is their potential of emphasizing, in the frequency domain, features that may be unnoticeable in the time domain. The oldest, and quite well known, is the Fourier transform introduced by Joseph Fourier [29].

This transform provides a decomposition of a signal into a basis of orthogonal trigonometric functions; due to the form of the Fourier basis functions, $e^{i\omega t}$, the spectrum is not localized in time and only the global frequency content of a signal is obtained [30]. This is acceptable regarding stationary and pseudo-stationary signals but quite unsatisfactory in the case of signals that are highly non-stationary, noisy, a-periodic etc. Note that, in this article, we focus on wavelets localized in time (and thus scales are time-scales by default), although space localization is also possible in general.

Under the circumstances, a methodology capable of obtaining the frequency content of a process locally in time is required [31]. An acceptable expedient to this condition of time–frequency localization was the Short-Time Fourier Transform (STFT) introduced by Gabor [32] in answer to his rhetorical question “…how are we to represent other signals, for instance a sine wave of finite duration?”.

The STFT comprises repeated FT calculations in a shifting, limited extent window, be it $g\left(t\right)$, of the signal centred at time t, which, in effect, partitions the time–frequency plane (also Phase Space, a loan from physics [31,33]) with rectangular cells of the same size and aspect ratio [30,34]. Thus, frequency and time information are both recovered from the aforementioned signal as each of the multiple FT spectra obtained corresponds to the time t.

The performance of the STFT analysis depends critically on the chosen window $g\left(t\right)$ since it determines the uncertainties in the frequency, $\Delta f$, and time, $\Delta t$, localization of the transform; furthermore, the product of these uncertainties is constant and of the order of unity. In the Appendix 9.3 of Gabor [32] the benefits of a Gaussian window $g\left(t\right)$ are shown as it minimizes the product of uncertainties ($\Delta f·\Delta t\ge 1/2$) to 1/2 (see also [35]); this type of STFT is called the Gabor transform. The caveat, in the case of the SFTF, is that once the function, $g\left(t\right)$, is selected, the uncertainties, $\Delta f$ and $\Delta t$, are both fixed, and the time and frequency resolution is constant in the time–frequency plane.

The disadvantages of the STFT, as discussed in the previous paragraph, stem from the uncertainty in frequency and time localization combined with the constant time and frequency resolution in the time–frequency plane. The circumvention of this problem needs a transform with basis functions, which may capture detail at different frequency scales at the same time.

One of the first modern methods that attempted to tackle this issue was proposed by Priestley [36] in 1965, under the frame of his Evolutionary Spectra methodology, by which he extended the classical Fourier transform with a generalized basis of orthogonal functions, indexed by both time and frequency. Even though this method was popular in the civil engineer community [37], it lacked a systematic way of producing these basis functions—something that severely impacted its usefulness. Another extended family of methods relies on Principal Component analysis or Singular Value decomposition techniques, and these are known as Empirical Orthogonal Function expansions (see [38] and the references therein).

These attempt to decompose a signal on an empirical basis of orthogonal eigenfunctions that are themselves derived from the signal itself, and thus there is no arbitrary choice of basis. Unfortunately these eigenfunctions not only have no frequency content but, in general, lack any intuitive physical meaning, and given the fact that the decomposition might not be unique, their interpretation becomes difficult. On the other hand, a method that is both easily applicable and readily interpretable is the wavelet transform, which also offers the major advantage of varying its time–frequency sensitivity, enabling dynamic adaptation to the desired time or frequency resolution.

The rest of this Section 2.1 is devoted to a brief outline of the evolution of the wavelet theory with emphasis on the continuous wavelet transform (CWT), the application to the study of geospace disturbances. According to Morlet et al. [39], the beginning is attributed to Gabor, and, in his works, the eponymous wavelet (see Section 2.2.1 Equation (2)) is mentioned as a Gabor wavelet. This transform employs almost the same base function as the Gabor STFT, which is a windowed complex exponential with a Gaussian window. In this case, the width of the window is adjustable, contrary to the STFT [40].

It would be compressed in time in order to obtain a higher frequency function or dilated to obtain a lower frequency function; these compressed or dilated functions would be shifted in time in order to obtain the required localization in frequency. The compressed wavelets are expected to provide high time resolution at some loss of frequency resolution, the dilated the opposite (see [22]). In brief, the wavelet transform overcomes the resolution limitations of the STFT by means of the constant relative bandwidth condition ($\Delta f/f$ = constant), [34]. We note, at this point, that, regardless of the similarity of the first wavelets to the Gabor STFT, many choices of the wavelet form are available depending on the context.

2.2. Wavelet Methods

2.2.1. Continuous Wavelet Transform (CWT)

The analysis of a function in time, be it $F\left(t\right)$, into an orthonormal basis of wavelets is conceptually similar to the Fourier Transform (FT); however, in the FT case, the spectrum is not localized in time as discussed in Section 2.1. A wavelet basis, on the other hand, is localized in both frequency f and time $\tau$, thus rendering itself suitable to the analysis and decomposition of non-stationary time series and transient signals [25,41,42]. The wavelet basis is derived by scaling and translating (shifting in time) one and only mother wavelet $\psi \left(t\right)$, to produce different daughter wavelets at each frequency and each temporal position.

The scaling is accomplished be means of a scale factor s, which transforms the mother function to $\psi (t/s)$, thus, changing its spread and hence its frequency. Even though it is tempting to interpret s as the inverse of the frequency ($f^{-1}$), it should be noted that the exact relation that maps s to the corresponding Fourier frequency is slightly more complicated and depends on the particular choice of the wavelet mother function, as well as the values of its parameters.

Due to this transformation, all daughter wavelets $\psi \left(\right(t-\tau )/s)$ are scale covariant and inherit from the mother wavelet zero average, square integrability ($\psi \left(t\right)$ in $L^2\left(R\right)$) and compact support or fast decay in order to warrant localization [23,30,43]. The wavelet transform of a function F(t) is calculated as the convolution of the function with the wavelet $\psi \left(\right(t-\tau )/s)$. The integral is replaced by a summation in the case $F\left(t\right)$ is represented by a discrete time series $F_i=F\left({\tau }_i\right)$; in both cases, the wavelet transform represents a mapping of $F\left(t\right)$ from the solely temporal to the time–frequency domain [44,45,46]:

$$W(\tau ,s)=\frac{1}{\sqrt{s}}·{\int }_{\infty }^{-\infty }dt·F\left(t\right)·{\psi }^{*}\left(\frac{\tau -t}{s}\right)=\sqrt{\frac{\delta t}{s}}·\sum _{i=1}^N{F}_i·{\psi }^{*}\left(\frac{{\tau }_i-t}{s}\right)$$

(1)

where * denotes complex conjugate and $\sqrt{s^{-1}}$ is necessary to satisfy the normalization condition that the wavelets must have the same energy at every time and scale (or frequency). In the discrete case, this normalization factor is extended to include the sampling time $\delta t$ of the series [41].

An analytic function $\psi \left(t\right)$ is classified as a wavelet if the following mathematical criteria are satisfied [40]:

1.

A wavelet must have finite energy, i.e., the integrated squared magnitude of $\psi \left(t\right)$ must be less than infinity.

2.

If $\Psi \left(t\right)$ is the Fourier transform of the wavelet $\psi \left(t\right)$, then $\Psi \left(t\right)$ must not have a zeroth frequency component ($\Psi \left(0\right)$ = 0), i.e., the mean of the wavelet $\psi \left(t\right)$ must equal zero. This condition is known as the admissibility constant.

3.

For complex wavelets, the Fourier transform $\Psi \left(t\right)$ must be both real and vanish for negative frequencies.

The WT allows for some freedom regarding the selection of the functional form of the mother wavelet, $\psi \left(t\right)$, yet restricted by the admissibility conditions mentioned above. Although a true physical time-series should be independent of the choice of the wavelets, it is highly recommended, for best results, to use a wavelet basis that resembles in form the signal [5,26]. As most astrophysical and geophysical time-series are usually composed of intermittent bursts of sinusoidal-like oscillations, it should not come as a surprise that the most common mother wavelet used is the Morlet wavelet [25,40,47,48], which consists of a complex plane wave modulated by a Gaussian (see Figure 1). The mother Morlet wavelet and the corresponding CWT of a discrete sequence $F_n$, see also Torrence and Compo [41] and Grinsted et al. [49], are:

$$\begin{array}{c}{\psi }_o\left(n\right)=\frac{1}{\sqrt[4]{\pi }}exp\left(i{\omega }_{o}n-\frac{n^2}{2}\right)\hfill \\ W_n^F\left(s\right)=\sqrt{\frac{\delta t}{s}}·{\displaystyle \sum _{n^{\prime }=0}^{N-1}}{F}_{n^{\prime }}{\psi }_o^{*}\left(\frac{n^{\prime }-n}{s}\delta t\right)\hfill \end{array}$$

(2)

where ${\omega }_0$ is the dimensionless frequency, and it is usually set to 6 to satisfy the admissibility condition [23], while $\eta$ is the dimensionless time [49]. The correspondence between the wavelet scale and Fourier frequency for Morlet is given by Equation (3), from which it is immediately derived that for ${\omega }_0=6$ the Fourier frequency is equal to $0.968/s$, and hence it could be argued that for this particular case $s\approx f^{-1}$ (Torrence and Compo [41], Table 1).

$$f=\frac{{\omega }_0+\sqrt{2+{\omega }_0^2}}{4\pi s}$$

(3)

The advantage of this wavelet being complex valued is that the CWT provides both amplitude and phase of the time series being analysed [50]. Furthermore, the similarity of the Morlet wavelet expansion to the Gabor transform, mentioned in Section 2.1, extends to the minimization of the product of uncertainties ($\Delta f·\Delta t=1/2$) [39]. An example of a wavelet spectrum utilizing the Morlet wavelet is shown in Figure 2.

The average of the wavelet power spectral density ${∥W_n^F\left(s\right)∥}^2$ is the global wavelet spectrum (panel c in Figure 2) and is given by:

$$\overline{W^F\left(s\right)}=\frac{1}{N}\sum _{n=0}^{N-1}{∥W_n^F\left(s\right)∥}^2$$

(4)

The global wavelet spectrum is an unbiased and consistent estimation of the true power spectrum of a time series and generally exhibits similar features and shape as the corresponding Fourier spectrum.

The need for tests that are able to distinguish statistically significant results from those due to random background fluctuations (significance tests) is an essential prerequisite to any interpretation attempt regarding the results obtained from the application of CWT (or FT) on a single data sample. Hence, an appropriate background model is required for this task; white and red-noise spectra are often used [41]. This is justified, as the time series of interest usually exhibit higher variability in the long wavelengths than in the short. Under the circumstances of the null hypothesis $H_o$ (purely random results), the time-series used in the wavelet spectrum calculation is generated by uni-variate lag-1 autoregressive–AR(1) or Markov–process.

$$F_n=\alpha ·F_{n-1}+z_n$$

(5)

where $\alpha$ weights the autoregressive term, $z_n$ is taken from Gaussian white noise (GWN) and power spectrum $P_k$ is [41,49]:

$$P_k=\frac{1-a^2}{1+{\alpha }^2-2\alpha ·cos\left(\frac{2\pi k}{N}\right)}$$

(6)

where k = 0 …N/2 is the frequency index.

As in [41,42] the test for $H_1$, significant at p or, equivalently within a confidence level of $1-p$, is performed based on the height of a wavelet power spectral peak above the background defined in Equation (5); the wavelet power spectral density ${∥W(t,f)∥}^2$ is shown in Torrence and Compo [41] to be ${\chi }_2^2$ distributed with two degrees of freedom due to the two parts, real and imaginary, of this WT transform. If the variance of the time series is equal to ${\sigma }^2$ the product ${\sigma }^2·P_F\left(f\right)$ is the theoretical wavelet power, which can then be compared with ${∥W(t,f)∥}^2$ for the computation of the confidence intervals for the peaks in a wavelet power spectrum (see [41,49]):

$$P\left(\frac{{∥W_n^F\left(s\right)∥}^2}{{\sigma }^2}\prec p\right)=\frac{{\chi }_2^2}{2}·P_k$$

(7)

where $P\left(X\right)$ indicates “probability of X”, and the significance is p where, after a time honoured tradition, $p=0.05$ for the 95% confidence interval.

In the example presented in panel (c) of Figure 2, the 95% significance level is depicted in the plot of the peaks of the global wavelet spectrum; likewise, in Figure 3. Lastly, it needs to be noted that the application of a ${\chi }_2^2$ based criterion presupposes a close to a Gaussian distribution of the time series $F_j$, [41]; in the case that the time series is far from Gaussian, some type of transformation must be considered, as suggested by Grinsted et al. [49].

The Continuous Wavelet Transform, presented in this section, is an effective method for the detection of the frequency and time structure of activity bursts, ubiquitous in most astrophysical and geophysical time-series. The study of the interactions of two or more time series, on the other hand, necessitates that this methodology is appropriately expanded; this is done in the following Section 2.2.2.

2.2.2. Cross Wavelet Transform (XWT) and Wavelet Coherence (WTC)

The Cross Wavelet Transform (XWT) extends the application of the CWT in the case of two or more time series (let them be X and Y); it underlines regions of coincident energy between signals and detects relative phase. It is, thus, a useful a wavelet–based method of investigation of the possibility of causal relationships in the time frequency space between them. If $W_n^X\left(s\right)$ and $W_n^Y\left(s\right)$ are the CWTs of X and Y, the cross-wavelet transform is defined as [40,41,42,49]:

$$W_n^{XY}\left(s\right)=W_n^X\left(s\right)·W_n^Y{\left(s\right)}^{*}$$

(8)

with $*$ denoting complex conjugate.

The result is, in general, complex, and its magnitude (or modulus) and phase (or argument) are examined separately; the modulus, $∥W_n^{XY}∥$, of the Wavelet Cross Spectrum indicates the covariance at scale $s(\approx f^{-1})$ and the non-dimensional time parameter n marks regions in the (t-f) space of high common power. The phase spectrum, obtained from $arg\left(W_n^{XY}\right)$, of the XWT represents the relative phase difference between the time-series to be compared [41,49]:

$${tan}^{-1}[\frac{Im\left(\left|W_n^{XY}\left(s\right)\right|\right)}{Re\left(\left|W_n^{XY}\left(s\right)\right|\right)}]$$

(9)

From the phase of the XWT the Wavelet Coherence, as a quantitative measure of the phase consistency (WTC) between $W_n^X$ and $W_n^Y$, will be derived in the remaining part of this section. It is noted in advance that the regions of high common power and consistent phase relationship indicate a possible causal relationship between X and Y.

The estimation of the significance of wavelet cross power peaks is distinctly similar to that of the wavelet spectral peaks outlined in Section 2.2.1. The red noise assumption, Equation (6), is used in the derivation of the null hypothesis and the estimates of the power spectra theoretical distributions $P_X\left(f\right)$ and $P_Y\left(f\right)$; the computation of the confidence intervals is provided in (Grinsted et al. [49] Equation (5)) (also (Torrence and Compo [41] Equation (31))):

$$P\left(\frac{∥{\mathrm{W}}_{\mathrm{n}}^{\mathrm{X}}\left(\mathrm{s}\right)·{\mathrm{W}}_{\mathrm{n}}^{\mathrm{Y}}{\left(\mathrm{s}\right)}^{*}∥}{{\sigma }_{{}^X}{\sigma }_Y}\prec p\right)=\frac{Z_{\nu }\left(p\right)}{\nu }·\sqrt{P_k^X{P}_k^Y}$$

(10)

where k is the frequency index and $\nu$ equals 1 if the wavelet in use is real valued and 2 if it is complex valued. The function $Z_{\nu }\left(p\right)$ is the confidence level associated with the significance p for a probability distribution function defined as the square root of the product of two ${\chi }^2$ distributions each corresponding to $P_k^X$ and $P_k^Y$ respectively; for example $Z_2(95\%)=3.999.$ and $Z_1(95\%)=2.182$ ((Torrence and Compo [41] Section 6c), see also [49]).

The Wavelet Coherence (WTC) is an estimator of the confidence level for each detection of a time–space region of high common power and consistent phase relationship between two time-series calculated by means of the Cross Wavelet Transform (Equation (8)). The measure of wavelet coherence is defined between two continuous wavelet transforms, and this may indicate coherence with high confidence level even though the common power is low; it closely resembles a localized correlation coefficient in time–frequency space and varies between 0 and 1. It is used alongside the Cross Wavelet Transform, as the latter appears to be unsuitable for significance testing the interrelation between two processes [52]. Following (Torrence and Webster [42] see Appendix) and Grinsted et al. [49], we define the wavelet coherence of two time series—let them be X and Y:

$$R_n^2\left(s\right)=\frac{{∥S\left({\mathrm{s}}^{-1}{\mathrm{W}}_{\mathrm{n}}^{\mathrm{XY}}\left(\mathrm{s}\right)\right)∥}^2}{S\left({\mathrm{s}}^{-1}{∥{\mathrm{W}}_{\mathrm{n}}^{\mathrm{X}}\left(\mathrm{s}\right)∥}^2\right)·S\left({\mathrm{s}}^{-1}{∥{\mathrm{W}}_{\mathrm{n}}^{\mathrm{Y}}\left(\mathrm{s}\right)∥}^2\right)}$$

(11)

where S is a smoothing operator and $R_n^2\left(s\right)\le 1$. The definition closely resembles that of a traditional correlation coefficient, which, in this case, is localized in the time frequency space. The smoothing operator S may be in the form:

$$S_W=S_{scale}\left(S_{time}{W}_n\left(s\right)\right)$$

(12)

where $S_{scale}$ denotes smoothing along the wavelet scale axis and $S_{time}$ smoothing in time. The smoothing operator is expected to have, by design, a footprint as small as the wavelet used (see (Grinsted et al. [49] Equation (10)) for the analytical expression of a smoothing operator for the Morlet wavelet example, which was initially described in Torrence and Webster [42]).

The statistical significance level of the wavelet coherence is estimated using Monte Carlo methods as described in Grinsted et al. [49].

3. Application of Wavelet Methods

The solar wind–magnetosphere coupling forces the geomagnetic field perturbations as it provides the necessary energy input, with the results manifesting themselves as geomagnetic activity. The variability in the solar magnetic field, interacting with the magnetosphere at multiple time scales, which correspond to different processes, is a major driver of this coupling. In this section, we discriminate the periodic behaviour of the geomagnetic field in large-scale, i.e., periodicities larger than 27 days (solar rotation) and short-scale, i.e., geomagnetic oscillations in the timescale of the order of seconds (Ultra Low Frequency waves–ULF).

3.1. Application of Wavelet Methods in Large–Scale Periodic Behaviour

3.1.1. Geomagnetic Field

As already mentioned in the introduction, the solar wind is modulated by the solar activity variations; its interaction with the terrestrial magnetosphere introduces the solar-terrestrial coupling, which appears as the geomagnetic activity dependence on solar activity fluctuations. Quantitative measures of the former are the geomagnetic indices, such as Ap, Kp, aa and Dst, as well as the auroral indices AE, AL and AU (see Mayaud [53] and Menvielle et al. [54] for a review). Many publications have been dedicated to the study of time-variations of geomagnetic activity indices.

Following the initial work of Bartels [55] discussing the time variation of geomagnetic activity indices Kp and Ap during 1932–1961, Fraser-Smith [56] used FFT analysis on a 38-year period (1932–1969) of Ap data and obtaining prominent peaks at 10.2, 7.04, 5.14, 4.10, 1.47, 1.09 and 0.5 years, as well as prominent peaks at 27.2, 27.6 days and moderate peaks at 54.0, 37.4, 30.5, 26.9, 18.7, 14.1, 13.7, 13.6 and 9.39 days.

As wavelet methods became popular in the 1990s, several articles were focused on the study of the interconnection of geomagnetic activity with the variation of solar energetic events and of solar wind parameters in various time-scales; the detection of shared periodicities provided evidence in support of this interconnection [10].

The ∼27-day periodicity, linked to the solar rotation, has been detected in all geomagnetic indices [10,57]. The second and first harmonic (∼9 and ∼13.5 days) have been also detected, and they were associated with the presence of two high speed streams per solar rotation at 1 AU [7].

Lou et al. [58] detected, using CWT, seasonal variations in the Ap index of 187, 273 and 364 days in the 1999–2003 period. They interpreted the annual and semi-annual quasi-periodic oscillations in terms of CME-magnetosphere interactions. For the explication of the 273 days (∼0.7 years) periodicity, on the other hand, they used estimates for periodic time-scales of equatorially trapped Rossby-type waves [3]. Using a significantly extended dataset, spanning four solar cycles in the 1966–2010 time period, Katsavrias et al. [10] detected time-localized common spectral peaks, between the fluctuations in the solar wind characteristics and the geomagnetic indices, which corroborated some of the previously found periodicities.

On the contrary, the prominent semi-annual periodicity detected in the Dst index (and less prominent in the Kp index) did not appear to have originated in the Solar wind parameter fluctuation (Figure 3). Therefore, it was attributed to reconnection processes in the day-side magnetosphere induced by the Russell and McPherron [59] effect. Furthermore, contrary to the Ap periodicity reported in Lou et al. [58], the AE index was found to have only an annual peak.

Nayar et al. [8] reported Quasi-Biennial Oscillations (QBO) of the geomagnetic activity index Ap; 1.3–1.4 years during even cycles and of 1.5–1.7 years during odd ones, likely propagating from the solar magnetic field (SMF) and were associated with double peaks of the solar cycle. These double peaks were attributed to a superposition of the usual 11-year cycle and wave trains with periodicities continuously varying from three years at solar maximum, to 1.7 years towards solar minimum [5].

A step further, Ou et al. [60] examined the second-order time derivatives derived from monthly means of the X, Y and Z field components recorded by the global network of ground-based observatories between 1985–2010. Having filtered out very high and very low frequency components (due to the use of the monthly mean and of the time-series derivative, respectively), they identified five principal periods at 1.3, 1.7, 2.2, 2.9 and 5.0 years.

The authors suggested that, even though QBOs in the geomagnetic field are generated from a common source (very high correlation between the geomagnetic QBO and the QBOs in solar wind speed and solar wind dynamic pressure), their features indicate that they primarily originate from the various current systems related to the solar wind-magnetosphere-ionosphere coupling process. The dependence of geomagnetic indices on the interplanetary magnetic field (IMF) polarities, toward and away (T and A, respectively) has been also investigated by utilizing wavelet analysis [61,62]. The results showed a significant N-S asymmetry of these indices, depending on the solar cycles.

In an attempt to link the quasi-periodic variations exhibited by the geomagnetic activity with the solar wind parameters, Andriyas and Andriyas [63] used wavelet analysis to study the interrelationship between ten solar wind coupling functions with the geomagnetic indices Dst and AL. The authors indicated that no single coupling function could explain the variances in geomagnetic indices at all the time-scales. Moreover, the authors observed that although there were periodicities common to the coupling functions and the AL index, the amount of correlation was poor when compared to the Dst index.

As already mentioned, the quasi-periodic behaviour of the geomagnetic activity is more complex than the that exhibited by solar activity parameters. This is due to the fact that the former is driven by two dominant drivers: (a) the high speed solar wind streams from coronal holes and (b) the intermittent solar activity (eruptive events). Feynman [11] was the first to compare aa index variations with sunspot numbers and, since the solar wind drives the geomagnetic activity, concluded that the 11-year solar cycle, as represented by the sunspot numbers, was very different from the one represented by the solar wind parameters and geomagnetic indices.

In particular, on annual scale, the geomagnetic index aa could be the result of two components: one originating from solar transient (or sporadic) activity and in phase with the solar cycle; the other was related to recurrent solar drivers with peak in the descending phase. The major driver of transient geomagnetic activity is the Interplanetary Coronal Mass Ejections (ICME). The solar recurrent activity, on the other hand, is driven by High Speed Solar Wind Streams (HSSWS) and Stream Interaction Regions (SIR).

Taking another step further, Katsavrias et al. [13] used cross-wavelet transform (XWT) and wavelet coherence (WTC) to study the relationship between transient and recurrent phenomena, i.e., ICMEs and CIRs and the corresponding magnetospheric response represented by geomagnetic indices Dst and AE, during the solar cycle 23. The authors verified that CIRs modulate the geomagnetic response during the ascending and descending phase, while ICMEs during the maximum of the cycle and the unusual active period of 2002–2005; nevertheless, this feature was evident in the 27-day periodic component only (Figure 4).

3.1.2. Magnetospheric Particles and Cosmic Rays

As already mentioned, the seasonal variation (annual and semi-annual) has long been recognized in geomagnetic activity. The discovery of the trapped particle populations in the inner magnetosphere (Van Allen radiation belts) opened a brand new scientific domain, which led to the launch of several scientific space missions in the 1990s. One of them, the Solar Anomalous and Magnetospheric Particle Explorer (SAMPEX) satellite, gave the appropriate measurements in order for the semi-annual variation (SAV) in the electron fluxes of the outer radiation belt to be discovered [64]. The nature of this variation were under debate for several years, as three possible mechanisms had been proposed:

1.

The axial effect [65], which is due to the variation of the position of the Earth in heliographic latitude.

2.

The equinoctial effect [66], as a consequence of the varying angle of the Earth’s dipole with respect to the Earth-Sun line.

3.

The Russell and McPherron [59] effect, a result of the larger z component of the interplanetary magnetic field (IMF) near the equinoxes in GSM coordinates, owing to the tilt of the dipole axis relative to the heliographic equatorial plane.

With the aid of wavelet methods Katsavrias et al. [67] linked the semi-annual variation to the Russell–McPherron effect: They used electron flux measurements from the Van Allen probes and Geostationary Operational Environmental Satellites (GOES) and performed cross-wavelet analysis with respect to the angles of the three proposed mechanisms. The calculated XWT and WTC spectral peaks indicated an in-phase relationship between the electron flux and the angle that controls the Russell–McPherron effect (Figure 5) at the ∼178 days periodicity. The same result was verified by other methods by Poblet et al. [68].

In addition, wavelet analysis revealed that the semi-annual variation was more pronounced during the descending phases of Solar cycles 22, 23 and 24 and, moreover, coexisted with periods of increased number of High Speed Streams (HSSWS) and decreased the Interplanetary Coronal Mass Ejections (ICME) occurrence, indicating that the SAV was a result of the modulation of reconnection produced by the variability of the controlling angles of the RM (and/or equinoctial) mechanism during periods of enhanced solar wind speed.

The aforementioned results were of significant importance to the outer radiation belt dynamics, since the angle, $\theta$, that controls the Russell–McPherron effect has been proven crucial for the prediction of the long-term variability of energetic electron fluxes [69].

Several large-scale periodicities have been also detected in the cosmic ray (CR) time-series. Galactic cosmic rays (GCRs) are energetic charged particles, spanning a broad range in energy (10${}^6$–10${}^{20}$ eV/nucleon), mostly with an out-of-solar-system origin. These high-energy GCRs can reach the Earth’s atmosphere and trigger particle and electromagnetic cascades, thereby, creating a great variety of secondary particles, such as neutrons, which can be measured by ground-based neutron monitors [51]. The transport of GCR in the heliosphere, as of all charged particles, is known to be influenced and controlled by the spatial and temporal distribution of the interplanetary magnetic field (IMF). Therefore, the GCR evolution is expected to share common periodic behaviour with solar/solar wind parameters.

Except for the well-known 11-year periodicity due to the solar cycle and the annual variations, peaks at 7.7, 5.5, 2 and 1.7 years have been found by Mavromichalaki et al. [70]. The authors used wavelet analysis to indicate that the period of 7.7 years was related to the 22-year cycle and consequently, with the polarity of the solar magnetic field. On the other hand, the period of 5.5 years was related to the enhanced power of the second harmonic, which arises from the asymmetric form of the 11-year solar cycle.

QBOs are also very prominent in neutron monitor data. In an extensive study using wavelet analysis, Kudela et al. [71] showed that, while a dominant quasi-periodicity centring at ∼1.7 years–was exhibited during the odd cycle 21 (1981-1984), there was an enhanced contribution from the 1.3 years quasi-periodicity during the even cycles 20 and 22.

Even though the aforementioned variations are linked with the solar activity dynamics, shorter-scale periodicities are linked with transient phenomena in the interplanetary space. Chowdhury et al. [72] detected significant periods in the range of 8–32 days including a solar rotational period of approximately 27 days, as well as prominent Rieger-type periodicities, during the maximum of solar cycle 24, associated with Solar eruptive events (Figure 6).

3.2. Application of Wavelet Methods in Short–Scale Periodic Behaviour

3.2.1. Ultra-Low Frequency Waves

Over the past decades, it has been well established that a geospace disturbance is the consequence of a geo-effective solar wind flow reaching the near-Earth space environment. These disturbances are associated with the variability of charged particles in the Earth’s radiation belts and the intensification of electric current systems with characteristic signatures on the geomagnetic field [73]. Periodic oscillations in the Earth’s magnetic field with frequencies in the range of a few millihertz (ultra low frequency ULF waves) can significantly influence radiation belt dynamics due to their potential for strong interactions with charged particle populations.

In particular, the ULF waves in the Pc4-5 band (1–25 mHz) can violate the third adiabatic invariant L${}^{*}$ of the energetic electrons. This drives radial diffusion by conserving the first two adiabatic invariants under the drift resonance condition $\omega =m{\omega }_d$, where $\omega$ is the wave frequency, m is the azimuthal wave mode number and ${\omega }_d$ is the electron drift frequency [74]. Radial diffusion is one of the most important mechanisms, since it has a dual role; it can lead to both energization [75,76,77,78,79] and loss of relativistic electrons [80,81,82].

ULF waves have been traditionally identified through the visual inspection of series of spectrograms based on the Fast Fourier Transform (FFT) or by the application of simple, automated or semi-automated (i.e., requiring some amount of human intervention) techniques [83,84]. Since the 1990s, wavelet spectral analysis has become popular, as it allows the quantitative monitoring of localized variations of power within the time series data [85]. Moreover, wavelet analysis can be superior to the Fourier spectral analysis when the spectral properties of transient, impulsive, short-lived or non-stationary signals need to be analysed, especially when the nature of the investigated signal is not a priori known.

An appropriate selection of the mother function allows for better detection of these signals and does so in a direct way, without the need for a moving window and the many arbitrary choices that arise from it, such as its length, step, envelope function and so on. This also produces outputs that are significantly smoother and thus allows for the better application of thresholds and more accurate identification of fast and short-scale phenomena. This is also evident by the series of studies that have highlighted the significance of applying wavelet analysis and especially its suitability for multipoint, small-scale disturbances, in the investigation of ULF wave events.

Wavelet analysis has been further used in the development of automated algorithms for the detection of ULF waves. Heilig et al. [86] developed an algorithm for the selection of possible upstream wave related pulsation events (frequencies in the 20–80 mHz range) from both ground and space magnetometer data. Balasis et al. [87] developed a time–frequency analysis tool for automated detection of ULF wave events in magnetic and electric field measurements made from multi-satellite missions and ground-based networks.

The Pc3 waves (in the 20–100 mHz frequency range) were identified by examining series of time–frequency spectrograms, produced with wavelet-based algorithms. This detection tool was later on implemented in a machine-learning-based model in order to detect ULF waves in the time series of the magnetic field measurements on board the low-Earth orbit CHAMP satellite [88]. Furthermore, Papadimitriou et al. [89] exploited this detection tool to create maps of Pc3-4 ULF wave power (in the 10–100 mHz frequency range–see also Figure 7) and a ULF wave index, which was then correlated with the solar wind parameters to elucidate the driving factors behind the onset and propagation of these waves.

The ability of the wavelet transform to characterize both the frequency and temporal content of waves makes it ideal for studying irregular pulsations (Pi1 and Pi2), as well. Especially, Pi2 ULF waves (7–25 mHz) are intimately connected with the onset of the expansion phase during magnetospheric substorms. Therefore, efforts have been devoted to the development of automated algorithms, utilizing a wavelet power threshold for Pi2 waves, for nowcasting/forecasting magnetospheric substorms from real-time ground-based magnetometer data [90] and ultraviolet images [91].

Wavelet analysis has also been applied in the detection of causal relations between the magnetosheath pressure and the waves observed in the magnetosphere. Archer et al. [92], investigated the impact of large amplitude, transient dynamic pressure pulses on the generation of compressional and poloidal mode waves (Pc5-6 frequency range) in the outer magnetosphere. They also used the cross-wavelet transform and wavelet coherence in order to detect the generation of well-defined standing waves (90${}^{\circ }$ phase relationship with coherence levels higher than 0.75). Using wavelet analysis, Archer et al. [93] showed that broadband jets in the subsolar magnetosheath directly and resonantly drive ULF waves in the magnetosphere at the so-called magic frequencies (roughly 0.7, 1.3, 1.9, 2.6, 3.3 and 4.8 mHz) as well as local field line resonances.

Another characteristic example of the use of wavelet methods in the magnetosheath–magnetosphere coupling is the work of Katsavrias et al. [94]. The latter authors used multipoint magnetic field observations to detect identical Pi2 pulsations in the magnetosheath (similar to the ones accompany substorm onsets), generated at the wake of a jet and in the outer magnetosphere. Using cross-wavelet analysis and wavelet coherence, showed that the phase difference between the two pulsations was in good agreement with the propagation time of a disturbance travelling with Alfv$\acute{e}$nic speed, indicating that they are directly related.

3.2.2. Radial Diffusion of the Trapped Electron Population in the Outer Radiation Belt

In the previous section, we discussed the importance of ULF waves, in the variability of the trapped electron population in the outer radiation belt, as drivers of the radial diffusion processes. The quantification of these processes is possible through the estimation of the radial diffusion coefficient (D${}_{LL}$), which represents the mean square change of the drift-shell (L${}^{*}$) for a large number of particles over time. According to Fei et al. [95], the D${}_{LL}$ is the sum of the effects of perturbations in the azimuthal electric field and the parallel magnetic field:

$$D_{LL}=D_{LL}^B+D_{LL}^E$$

(13)

These two components of the radial diffusion coefficients are given by:

$$D_{LL}^B=\frac{{\mu }^2{L}^4}{8q^2{\gamma }^2{B}_E^2{R}_E^4}·\sum _m^{}{m}^2{P}_m^B\left(m{\omega }_d\right)$$

(14)

$$D_{LL}^E=\frac{L^6}{8B_E^2{R}_E^2}·\sum _m^{}{P}_m^E\left(m{\omega }_d\right)$$

(15)

where $\mu$ is the first adiabatic invariant, L is the Roederers L${}^{*}$, q is the charge of the diffused electrons, $\gamma$ is the Lorentz factor, $R_E$ is the Earth’s radius and $B_E$ is the strength of the equatorial geomagnetic field on the Earth’s surface. Moreover, P corresponds to the wave power at a specific drift frequency (${\omega }_d$) for all the azimuthal mode numbers (m).

It is clear, from this formulation (Equations (14) and (15)), that accurate calculations of the D${}_{LL}$ require equally accurate spectral analysis of the magnetic and electric field measurements used in the computation of the Power Spectral Density (PSD) of the Pc4-5 ULF waves. Since the temporal evolution of the PSD is of outmost importance for the calculation of the D${}_{LL}$, wavelet analysis is one of the most commonly used methods [96,97,98,99,100].

4. Revisiting the Estimation of Wavelet PSD and Comparison with FFT

One of the most often-heard criticisms of wavelet analysis is that, for most cases, the mother wavelets do not form an orthonormal basis, meaning that their inner product for frequencies at different octaves, i.e., at frequencies with ratios that are integer powers of two, is not zero [101]. Even though, this is technically true, for a careful choice of the mother function and its parameters, the derived wavelets can be quasi-orthogonal, in the sense that even if their product is not exactly zero, it is a very small number that can be safely approximated by zero.

As an example, by choosing the value of 6 for the ${\omega }_0$ parameter of the Morlet wavelet (Equation (2)), the inner product of the generated wavelets is smaller than ${10}^{-6}$ and so can be considered zero for all practical applications.

A second argument against wavelets is that the derived wave power depends on the values of the parameters employed in the methodology—something that makes comparisons against other methods and particularly against spectra derived by means of the FFT problematic. This becomes even more important when the results of the analysis are being used as inputs for the calculation of other quantities, such as the radial diffusion coefficients discussed in Section 3.2.2 above. This is a critical issue as the differences between even small changes in the parameters can be quite substantial; however, it can also be alleviated by applying the proper normalization technique. Note that the corresponding code for the estimation of PSD using this normalization technique is publicly available at https://synergasia.uoa.gr/modules/document/index.php?course=PHYS120 accessed on 9 March 2022.

To examine all these and better illustrate this method, we produce an artificial signal, to be used as the benchmark for the tests to follow. The signal is a modulated sinusoidal, with a time-varying frequency and some additive noise and is given by:

$$x\left(n\right)=10·exp({-\frac{(t\left(n\right)-t_0)}{N-1}}^2)·sin\left(2\pi t\right(n)/f(n\left)\right)+R\left(n\right)$$

(16)

with $f\left(n\right)$ varying linearly from 1 up to 125 mHz, $R\left(n\right)$ being Gaussian random noise with a mean of zero and standard deviation of 1 and the exponential term signifying a Gaussian envelope in the amplitude of the sinusoidal. The signal was constructed as a series of 20,000 points with a sampling time of $\delta t$ = 4 s, so that $t\left(n\right)=1+(n-1)·\delta t$ with n from 1 up to 20,000 and can be seen at the top panel of Figure 8.

To compute the wavelet spectrum of this signal, the Torrence & Compo code was used, with the Morlet mother function and an ${\omega }_0=6$. To define the scales (and correspondingly the frequencies) on which the analysis will be performed, we used the convention in the same publication [41] to construct 37 scales, starting from $s_0=2·\delta t$, which practically corresponds to the Nyquist frequency for this test signal and using a scale step dj of 0.25, implying that we generate four scales for each octave. The scales themselves follow the relation $s_m=s_0·2^{m·dj}$, with m from 0 up to 36, which yields scales from $2·\delta t=8$ sec up to $2^{12}·\delta t=$16,384 s. The wavelet power $W_n\left(f\right)$ of the signal under these parameters is shown in the middle panel of Figure 8, from which we can see both the rising tone characteristic of this signal and its initially increasing and afterwards diminishing amplitude.

To compare against Fourier spectra, we perform a windowed-FFT analysis, by using a moving window of 900 points, i.e., one hour in duration, moving it by a step of only one point each time and computing the Fourier spectrum on the same range of frequencies as for the wavelet case, taking into account the relation between the two as given in Equation (3). The orthogonality of the Fourier transform and the conservation of the signal energy that it implies, allows us to compute the sum of the squares of the magnitude of the Fourier coefficients, along this frequency range and obtain a total wave power for each point in time. Doing naively the same for the wavelet case produces results that differ significantly from the FFT ones.

As the frequencies themselves can be chosen almost arbitrarily for the wavelet case, we can obtain different results even using a different frequency resolution dj, since the summation will now extend to more frequencies on which the wavelet representation is highly correlated. All these points are easily illustrated in the bottom panel of Figure 8, where the blue and black curves correspond to the sum of the wave power for the wavelet method (one with dj = 0.25 and the other of dj = 0.5), while the red curve to the result of the windowed FFT analysis. Thankfully, the normalization that is given in (Torrence and Compo [41] Equation (14)) provides us with the way to re-introduce the conservation of the signal’s energy in the wavelet formalism, so now the Total Wave Power for each moment in time can be computed as:

$$T_n=\frac{dj\delta t}{C_{\delta }}\sum _j\frac{|W_n\left(s_j\right){|}^2}{s_j}$$

(17)

with $C_{\delta }$ an empirically derived factor that for the Morlet case (with ${\omega }_0=6$) is equal to 0.776. By following this normalization, the total power along all frequencies examined in this example produces the green curve, which matches perfectly the red curve, which was the result of the windowed FFT analysis. Hence, by revisiting and adopting a normalization approach suggested more than 20 years ago (but overlooked in many modern wavelet studies involving space physics data), we are able to conserve the energy of the time series under analysis in the time–frequency domain using the wavelet transform.

Applying this normalization allows the full arsenal of the wavelet methodology and all its benefits to be used in a robust and consistent way for the estimation of the wavelet power spectral density within a specific frequency range and with respect to time, so that other parameters that are based on it, e.g., the radial diffusion coefficients, can now be accurately derived and compared with those that are procured from more traditional, Fourier-based methods. Especially for the radial diffusion coefficients case, the calculation of the wavelet power spectral density that is needed for their calculations can now immediately be achieved by dividing the total wave power by the step in frequency range $\delta f$, hence $P_n=T_n/\delta f$.

5. Summary

In this review, we highlighted the significant contributions of the wavelet methods in the research of solar-terrestrial coupling and geospace disturbances. Wavelet analysis can be superior to other methods for spectral analysis since its time-localization allows for the better detection of non-stationary signals (especially when their nature is not a priori known), without the need for a moving window and the many arbitrary choices that arise from it.

Since, the 1990s, wavelet methods have been applied in the accurate identification of not only fast and short-scale phenomena, such as ULF waves (geomagnetic oscillations in the timescale of the order of seconds) and the subsequent radial diffusion coefficients, but also in the identification of a variety of large-scale phenomena (periodicities larger than 27-day solar rotation) in the Sun, the solar wind, the geomagnetic field, the trapped particle population in the inner magnetosphere and cosmic rays.

Finally, we presented an optimal approach for the accurate estimation of the PSD that is compatible with the results obtained from FFT analysis, highlighting how wavelets can be used consistently regardless of the choice of mother function or other parameters. The corresponding code for the estimation of PSD using the aforementioned method is publicly available at https://synergasia.uoa.gr/modules/document/index.php?course=PHYS120 accessed on 9 March 2022.