Influence of Solar Wind on Secondary Cosmic Rays and Atmospheric Electricity

A relationship between the heliospheric magnetic field, atmospheric electric field, lightning activity, and secondary cosmic rays measured on the high mount of Lomnický Štít (2,634 m a.s.l.), Slovakia, during the declining phase of the solar cycle 24 is investigated with a focus on variations related to solar rotation (about 27 days). The secondary cosmic rays are detected using a neutron monitor and the detector system SEVAN, which distinguishes between different particles and energies. Using spectral analysis, we found distinct ∼27-day periodicities in variations of B_x and B_y components of the heliospheric magnetic field and in pressure-corrected measurements of secondary cosmic rays. The 27-day variations of secondary cosmic rays, on average, advanced and lagged the variations of B_x and B_y components by about 40° and −140°, respectively. Distinct 27-day periodicities were found both in the neutron monitor and the SEVAN upper and middle detector measurements. A nondominant periodicity of ∼27 days was also found for lightning activity. A cross-spectral analysis between fluctuation of the lightning activity and fluctuation of the heliospheric magnetic field (HMF) showed that fluctuation of the lightning activity was in phase and in antiphase with B_x and B_y components of the HMF, respectively, which is in agreement with previous studies investigating the influence of solar activity on lightning. On the other hand, the ∼27-day periodicity was not significant in the atmospheric electric field measured in Slovakia and Czechia. Therefore, no substantial influence of B_x and B_y on the atmospheric electric field was observed at these middle-latitude stations.

Introduction

It is well known that the heliospheric magnetic field (HMF) influences the intensity of cosmic rays (CRs) entering the Earth’s atmosphere (Usoskin et al., 1998; Owens and Forsyth, 2013). Primary (galactic) cosmic rays consist of high-energy particles, mainly protons with energies of ∼GeV and higher that interact with the Earth’s atmosphere, creating a shower of different species of charged particles, neutrons, and gamma ray photons. These showers, also called secondary cosmic rays (SCRs), are usually measured using neutron monitors (NMs). The count rates measured using NMs depend not only on the intensity of CRs but also on the strength and orientation of the geomagnetic field and the altitude of the measuring site. Therefore, many NMs are located at high latitudes or high altitudes where the so-called cutoff rigidity is lower, which means that larger count rates are measured as particles with lower momentum can reach such sites (Shea and Smart, 2000; López-Comazzi and Blanco, 2020). It was found that the intensity of CRs is anticorrelated with solar activity, exhibiting about 11-year cycle variations (Usoskin et al., 1998). The lower the solar activity, the higher the intensity of CRs that is observed because more galactic CRs penetrate the heliosphere and reach the Earth. Voiculescu and Usoskin (2012) and Voiculescu et al. (2013) showed using statistical studies that solar activity and the HMF might impact cloud cover in specific regions on long timescales. Variations of CRs in the Earth’s atmosphere are important from many aspects which are as follows: a possible link between climate and intensity of CRs via potential influence of CRs on cloud condensation nuclei (Kristjánsson et al., 2008; Kirkby 2008; Svensmark et al., 2009; Gray et al., 2010), the dominant role of CRs in the production rate of isotope ¹⁴C used for dating (e.g., Gosse and Klein, 2015), modulation of doses received by humans and electronics on board aircrafts (Yang and Sheu, 2020), and possible influence on seed electrons for lightning initiation (Dwyer and Uman, 2014).

Periodicities or quasiperiodicities shorter than the 11-year period of the solar cycle were also identified in the CR intensity. For example, Kudela et al. (2002), Kudela et al. (2010), and Chowdhury and Kudela (2018) focused on periods longer than the solar rotation and studied Rieger-type fluctuations around ∼154 days and quasiperiodicities around ∼1.3, 1.7, and 2.25 years. The intensity of these periodicities usually varied during the solar cycle, and their origin is still under debate (Bazilevskaya et al., 2014). The Rieger-type fluctuations (∼154 days) mainly occur around the solar maximum and were first investigated by Rieger et al. (1984). On the other hand, Saad Farid (2019) and López-Comazzi and Blanco (2020) also analyzed shorter periods such as the period of solar rotation (∼27 days) and its harmonics (mainly ∼13.5 days) besides the quasiperiodicities around 130 days and the Rieger-type fluctuations and quasiperiodicities in the range of 250–470 days.

The period of solar rotation (∼27 days), which is the main focus of this study, arises from the tilt of the solar magnetic dipole with respect to the rotation axis, coronal-hole high-speed streams, and corotating interaction regions characterized by the interaction of fast and slow solar wind (Gosling and Pizzo, 1999; Grieder, 2001; López-Comazzi and Blanco, 2020). The much weaker fluctuations at the period of ∼13.5 days might be just the second harmonics; however, the presence of the heliospheric current sheet (HCS) that separates oppositely oriented field lines of the HMF and occurs two times during the solar rotation (Owens and Forsyth, 2013) might also contribute to the observed fluctuations of SCRs at ∼13.5 days. In an ideal Parker spiral (Parker, 1958), the polarity of the dawn–dusk component B_y and the sunward component B_x of the HMF changes from away from the Sun (B_x < 0 and B_y > 0) to toward the Sun (B_x > 0 and B_y < 0), exhibiting one cycle during one solar rotation; the geocentric solar ecliptic (GSE) coordinate system is considered here. Interestingly, a dependence of the atmospheric pressure and atmospheric electric field on the polarity of B_y was found in the polar regions, mainly in the winter Antarctic when the air was relatively isolated from the rest of the world (Mansurov et al., 1974; Burns et al., 2008; Lam and Tinsley, 2016). It is assumed that the pressure dependence on B_y results from an external electric field that is superposed at high latitudes on the internal electric field, which is maintained by thunderclouds in the atmospheric global electric circuit (Rycroft et al., 2000). The external electric field occurs because of the relative motion and the interaction of the solar wind (SW) with the Earth’s magnetosphere, and its value in the SW sensed in the Earth frame (GSE coordinates) is given by the vector product E = -v x B, where v is the speed of the SW and B is the HMF, respectively. The component of the external electric field parallel with the Earth’s axis is then $E_z=-\left(v_x{B}_y-v_y{B}_x\right)\approx -v_x{B}_y$ as $\left|v_x\right|\mathrm{>}\mathrm{>}\left|v_y\right|$. The relative fluctuations of B_y are much larger than the relative fluctuations of v_x; consequently, E_z is mainly controlled by B_y (Lam and Tinsley, 2016). The E_z component is directed along the Earth’s axis and does not change during the Earth’s rotation if B_y is constant (Lam and Tinsley, 2016). A dominant component of the external electric field perpendicular to the Earth’s axis is the dawn–dusk component $E_y=-\left(v_z{B}_x-v_x{B}_z\right)\approx v_x{B}_z$ which exhibits 1-day periodicity due to the Earth’s rotation. The E_z component, controlled by B_y, dominates in the polar regions. It maps along the Earth’s magnetic field lines to the ionosphere and contributes to the total potential difference between the ionosphere and the Earth at high latitudes (Burns et al., 2008). The mechanism by which the atmospheric pressure is influenced by this external electric field is, however, unknown. It is hypothesized that an increase in the atmospheric vertical current that is associated with the increase in the potential difference between the ionosphere and the Earth’s surface influences the cloud microphysics, especially at cloud boundaries. This might result in different growths of condensation nuclei and hence in different radiation properties (opaqueness) of clouds or release of latent heat (Nicoll and Harrison, 2009; Gray et al., 2010; Lam and Tinsley, 2016). No direct dependence of atmospheric pressure and the electric field on B_y has been observed at middle or low latitudes. Variations of the electric field associated with global thunderstorm activity, local weather, cloud electrification, and changes of air conductivity owing to fluctuations of aerosol concentration dominate here. However, several studies have discussed how lightning occurrence over Great Britain and Japan likely depends on the B_y polarity and the solar rotation period (Owens et al., 2014; Owens et al., 2015; Miyahara et al., 2018).

This study presents measurements of SCRs taken using an NM and the Space Environment Viewing and Analysis Network (SEVAN) located at the mountain peak of Lomnický Štít (LS) in the High Tatras, Slovakia, and their comparison with SW and HMF data and nearby measurements of the atmospheric electric field. The aim is to search for similarities in these data series, especially similarities related to the period of solar rotation, and to discuss a potential influence of the HMF and SCRs on atmospheric electricity.

Measurement Setup

The near-Earth components of the SW velocity (v_x, v_y, and v_z) and the B_x, B_y, and B_z components of the HMF used in this study were downloaded from NASA/GSFC’s Space Physics Data Facility’s OMNIWeb service (https://omniweb.gsfc.nasa.gov/).

The SCRs are measured using two different types of detectors located on the peak of LS (49.195°N, 20.213°E) at the altitude of 2,634 m. The first is an NM installed on the top of LS in 1981. It is of the NM-64 type and consists of eight SNM-15 counting tubes with a length of 200 cm and a diameter of 15 cm, filled with BF₃ (Kudela and Langer, 2009). The NM at LS is a part of the global neutron monitor, which means that the data are included in the Neutron Monitor Database (NMDB). Cutoff rigidity of the NM at LS is 3.84 GV. The second is the SEVAN system, which has three separated plastic scintillators as detectors (Chilingarian et al., 2018). The upper and lower scintillators have the dimensions 100 × 100 × 5 cm and are separated from the middle scintillator by 4.5-cm-thick and 100 × 100 cm–wide lead absorbers. The middle detector is composed of five plastic slabs and has the dimensions 50 × 50 × 25 cm. The upper detector—channel 1—is mainly sensitive to low-energy charged particles, mainly electrons, and partly also to gamma rays (Chilingarian et al., 2018). The estimated energy threshold, considering the roof above the detector, is 7–8 MeV (Kudela et al., 2017). The middle scintillator—channel 2—mainly detects neutrons or high-energy gamma photons. It might be partly sensitive to electrons with energies of several tens of MeV. The lower detector—channel 3—registers high-energy charged particles, mainly muons with energies exceeding ∼250 MeV. Besides the count rates in individual channels, different combinations of all three channels are also evaluated. Specific combinations can be used to identify the type of incident particles. For example, combinations [1,0,0], [0,1,0], and [1,1,1] mean that a particle was registered only in channel 1, only in channel 2, and in all channels simultaneously, respectively. The combination [1,0,0] is mainly used for detection of electrons or positrons (partly also photons) with energies around ∼10 MeV. On the other hand, the combination [0,1,0] provides information about detection of neutrons, and the [1,1,1] combination provides information about the detection of high-energy (>250 MeV) charged particles, mainly muons (Chilingarian and Reymers, 2008; Chilingarian et al., 2021).

The electric field is measured on LS using an electric field mill EFM 100 sensor made by the Boltek company (Kudela et al., 2017). Relatively large values of the electric field are measured on the top of LS as it is a sharp, rocky mountain peak (Kudela et al., 2017; Chum et al., 2020). It should be noted that the mountain of LS represents a relatively conductive material, relative to the ambient air, and therefore locally enhances the intensity of the electric field. The EFM 100 at LS is in operation only during summer to prevent damage to its rotor due to the frequent ice coverage in winter. Therefore, the data of the electric field measured in Czechia at the stations of Panská Ves (50.527°N, 14.568°E, altitude 318 m) and Studenec (50.258°N, 12.518°E, altitude 666 m) are also provided for completeness and comparison. These stations are the closest available stations providing electric field data at a similar latitude and are part of a global atmospheric electricity monitoring network (Nicoll et al., 2019). More precisely, the electric field data are given as potential gradient (PG), which means that positive values correspond to the downward pointing electric field (electrons are accelerated upward). All the EFM 100 sensors used are installed in inverted positions to minimize precipitation noise.

To evaluate lightning activity, we used data provided by the World Wide Lightning Location Network (WWLLN), which operates about 70 sensors in the frequency range from 3 to 30 kHz (Rodger at al., 2004). Localization of individual discharges is based on the detection of sferics, impulsive signals radiated by lightning and propagating in a waveguide formed by the Earth’s surface and the bottom of the ionosphere. The WWLLN provides the exact times and locations of the detected lightning together with their estimated energy, energy error, and number of the stations used for the estimation of the relevant locations and energies. In this study, we used lightning detections which occurred between latitudes 48° and 51°N and longitudes 11° and 23°E. This region roughly covers the area of Czechia and Slovakia. The detection efficiency of the WWLLN did not change during the analyzed period.

Methods

To investigate the influence of the SW and the HMF on SCRs and the atmospheric electric field, we search for similar fluctuations in the HMF, SW, SCRs, and atmospheric electric field by performing spectral analysis, namely, Fourier transform and wavelet transform (WT). Because various quantities have different absolute values and units, it is useful to work with normalized unitless quantities. All the quantities are therefore normalized by applying Eq 1 before performing the spectral analysis as follows:

$a_{norm}=\frac{a-mean\left(a\right)}{{\sigma }_a},(1)$

where a is the analyzed quantity (components of SW velocity, components of the HMF, counts measured by the NM and the SEVAN, and the PG), and σ_a is the standard deviation of the distribution of the quantity a. The normalization allows for an easy comparison of fluctuations and the spectral content of various quantities.

All the data are available with 1-min resolution. It is necessary to consider that the data might contain outliers. For example, enhancements of SCRs related to the large atmospheric electric field were observed during thunderstorms (Kudela et al., 2017; Chum et al., 2020). The extreme values are therefore first removed from the 1-min data before the application of Eq 1. The extreme values or outliers were removed following the approach taken by López-Comazzi and Blanco (2020). It means that all the values lying outside the range [Q₁-1.5(Q₃-Q₁), Q₃+1.5(Q₃-Q₁)] were considered outliers, where Q₁ and Q₃ are the first quartile and the third quartile, respectively. The extreme values were removed for all the data, except the lightning counts in the considered area (48–51°N, 11–23°E). To reduce data volume and to speed up the data processing, we computed 2-h mean values from the original 1-min data. It should be noted that we are interested in fluctuations in the order of days. As the 2-h mean values act as a low-pass filter with a cutoff period of 2 h, the 2-h mean values were calculated with a 1-h step to fulfill the sampling theorem (sampling frequency should be at least twice the highest frequency in the signal). The mean values were only calculated if more than 50% of the data were available in the given 2-h intervals. Otherwise, the values were substituted by linear interpolation before the spectral analysis.

It should also be noted that only a rough proxy for fair-weather electric field data is obtained by removing extreme values from PG data using the method mentioned in the previous paragraph. The PG data might be influenced by the presence of moderately charged clouds, fog, or aerosols that significantly change the air conductivity and hence the PG. However, this approach can be justified as we are interested in the spectral content of the PG signal, rather than in absolute values. It is necessary to have as few data gaps as possible for spectral analysis. The main aim is to search for periodicities related to solar rotation, and it is unlikely that periodicities of fog or aerosol concentrations would follow the periodicities of solar rotation. Moreover, we found that we obtained a similar Carnegie curve for the Panská Ves station to that obtained in the study by Nicoll et al. (2019), who used another, stricter proxy for fair-weather data.

We chose the period from 2016 to 2019 for the study as for this period, all the data of interest are available for most of the time.

The measured values of SCRs depend on atmospheric pressure, as is documented in Figure 1A–D, which shows scatterplots between the counts recorded in SEVAN channel 1 ([1,0,0] combination), channel 2 ([0,1,0] combination), combination [1,1,1], and in the NM. To remove the pressure influence on the measured values of SCRs, the counts were first corrected on pressure using the linear relationships displayed in the individual plots in Figure 1 by green lines. Actually, the relationship between the counts and atmospheric pressure is expected to be exponential. However, the scatterplots indicate that linear correction does not introduce a substantial error in the observed range of pressure variations. Obviously, the counts for combination [1,1,1] (Figure 1C) exhibit the most complicated fluctuations and pressure correction is the least effective. The residual counts, obtained by subtracting the linear regression curves (green lines) from the measured counts, are used for spectral analysis after their normalization using Eq 1.

FIGURE 1

FIGURE 1. Scatterplots between measured counts (2-h mean values of 1-min data) and atmospheric pressure on LS. Green lines indicate linear regressions. (A) SEVAN [1,0,0] combination, (B) SEVAN [0,1,0] combination, (C) SEVAN [1,1,1] combination, and (D) NM.

In addition to Fourier and wavelet transform, cross-spectral analysis is performed for the quantities that exhibit similar spectra, namely, a distinct peak around the period of solar rotation. The cross-spectra make it possible to investigate a phase relation between quantities of interest. For example, cross-spectral densities (CSDs) for normalized pressure-corrected NM counts N_NM and the normalized B_x component of the HMF are calculated using relation Eq. 2 as follows:

$CSD\left(\omega \right)=N_{NM}\left(\omega \right)\cdot B_x^{\ast }\left(\omega \right)=\left|N_{NM}\left(\omega \right)\right|\cdot \left|B_x\left(\omega \right)\right|\cdot e^{i\left[{\phi }_{NM}\left(\omega \right)-{\phi }_{Bx}\left(\omega \right)\right]},(2)$

where $N_{NM}\left(\omega \right)$ is the Fourier transform of the normalized residual NM counts, $B_x^{\ast }\left(\omega \right)$ is the complex conjugate of the Fourier transform of the normalized B_x component of the HMF, and i² = -1. The phase shift $\varphi \left(\omega \right)$ between the two signals at a given frequency is then obtained as follows:

$\varphi \left(\omega \right)={\phi }_{NM}\left(\omega \right)-{\phi }_{Bx}\left(\omega \right).(3)$

The NM counts advance the B_x variations for positive $\varphi \left(\omega \right)$ and vice versa.

Results

Heliospheric Parameters and Cosmic Rays

Figure 2A,B shows the components of the HMF, together with its absolute value B, and the components of SW velocity (Figure 2C,D), including its absolute value v, for the years 2016–2019 in GSE coordinates. The B_x and B_y components have similar amplitudes and are mostly anticorrelated, which is an expected situation for the Parker spiral (Owens and Forsyth, 2013). The dominant component of the SW velocity is the v_x component (Figures 2C, D). The mean value of v_x in the GSE coordinates calculated over the displayed period of time is –425 km/s. The negative sign means that the SW velocity is directed from the Sun to the Earth. Figure 3 shows the power spectral densities (PSDs) of the normalized components of the HMF and SW velocity (normalized quantities of those displayed in Figure 2). The PSDs are displayed as a function of the period, rather than of the frequency. Only periods up to 40 days are presented. Figure 3A shows that the PSDs of B_x and B_y have a distinct peak around the period of 27.5 days. Both peaks have similar values, and the peak for B_x (red) is only slightly larger than that for B_y (blue). Figure 3B demonstrates that no significant peak around the period of solar rotation (∼27 days) was observed for the B_z component. Figure 3C shows that the PSD of the v_x component of SW velocity has two distinct peaks associated with solar rotation, namely, at the periods of 25.7 and 27 days. A likely explanation for these two peaks is a differential rotation of the Sun (Owens and Forsyth, 2013). SW originating from coronal holes around the solar pole exhibits a period of around 25 days, whereas SW arriving from the equatorial regions of the Sun is characterized by a period of around 27 days. Obviously, the periodicity of v_x propagates into the periodicity of the absolute value of SW velocity v (Figure 3D).

FIGURE 2

FIGURE 2. Components of the HMF (A, B) and components of SW velocity (C, D) in GSE coordinates.

FIGURE 3

FIGURE 3. Power spectral densities (PSD) of normalized components of the HMF (A, B) and velocity of SW (C, D).

Figure 4 shows time fluctuations of the normalized pressure-corrected counts recorded in SEVAN channel 1 ([1,0,0] combination), channel 2 ([0,1,0] combination), combination [1,1,1], and in the NM. To ensure homogeneity of the data, only counts from June 10, 2016 to April 15, 2019 are used for the analysis. Before this selected interval, changes in the acquisition electronics of the SEVAN were made. After this interval, many data gaps and discontinuities in SEVAN data occurred. Figure 5 presents the PSD of the normalized counts that are displayed in Figure 4. The PSD of SEVAN [1,0,0] and [0,1,0] combinations and the PSD of NM counts exhibit remarkable peaks at the period of solar rotation, around 27 days. However, no such distinct peak was observed for the SEVAN [1,1,1] combination (Figure 5C) that detects high-energy charged particles (muons). A likely explanation for that is that the HMF and SW are not able to sufficiently modulate the CR particles of high energies that are responsible for muon production because of the large gyro-radii of high-energy particles. There is a relatively good similarity between the spectra obtained for the B_x and B_y components of the HMF and the PSD obtained for the SEVAN [0,1,0] combination and NM counts. Both these detectors are mainly sensitive to neutral particles, the SEVAN [0,1,0] combination being also partly sensitive to high-energy gamma rays. A partial similarity with the PSD of B_x and B_y is also observed for the SEVAN [1,0,0] combination, sensitive to low-energy charged particles (approximately above 7 MeV), electrons and positrons. A partial detection of gamma rays cannot be excluded.

FIGURE 4

FIGURE 4. Normalized pressure-corrected counts observed using the SEVAN (A, B, C) and an NM (D) on LS.

FIGURE 5

FIGURE 5. Power spectral densities (PSD) of normalized pressure-corrected counts measured using the SEVAN (A, B, C) and an NM (D) on LS.

It is useful to investigate the time evolution of the spectra. Figure 6A–F presents the wavelet transform of the normalized B_x, v_x, and normalized counts recorded using the SEVAN [1,0,0], [0,1,0], and [1,1,1] combinations and the NM, respectively. The peak around the period of ∼27 days in the WT of B_x is most distinct in 2016 and in the beginning of 2017, which corresponds to the declining phase of the solar cycle 24. A second harmonic, ∼13.5 days, is also clearly visible in the WT of B_x. The time evolution of the second harmonic, however, does not follow the intensity observed at the period of ∼27 days. An almost identical result was obtained for the WT of B_y (not shown). Periods of ∼27 and ∼13.5 days can also be clearly identified in the WT of normalized v_x. Intensities observed at the periods of ∼27 and ∼13.5 days are comparable, which can also be seen in the PSD displayed in Figure 3C. It should be noted that a frequency resolution of the WT is not sufficient to distinguish the periods of 25.7 and ∼27 days observed in Figure 3C. Figures 6C, D, F show that the fluctuations of SCRs measured using the SEVAN [1,0,0] and [0,1,0] combinations and the NM exhibit relatively broad peaks around the period of ∼27 days. These peaks are more or less randomly distributed during the analyzed period of time. Intervals with relatively high intensities (distinct peaks) are followed by intervals with lower intensities. Therefore, it is likely that the fluctuations of SCRs are also modulated by other mechanisms besides the HMF and SW velocity. The exact origin of these mechanisms is unknown. We note that the counts were corrected to the local ambient pressure. However, we cannot exclude that a pressure or temperature profile in the atmosphere above the SEVAN and the NM could modulate the recorded counts (Riádigos et al., 2020). A potential influence of some other, unconsidered parameters of the solar and geomagnetic activity and of the atmospheric state can also not be excluded.

FIGURE 6

FIGURE 6. Wavelet transform of (A)B_x component of the HMF, (B)v_x component of SW velocity, and (C–F) normalized pressure-corrected counts measured using SEVAN combinations [1,0,0], [0,1,0], and [1,1,1] and the NM, respectively, on LS.

Atmospheric Electric Field

Figure 7 presents the normalized PG recorded at Studenec, Panská Ves, and LS. As mentioned in the Introduction, measurements on LS were not performed in winter to prevent damage from icing on the rotor of the EFM 100 sensor. It can be noticed that a higher PG was usually observed in winter than in summer. Large concentrations of aerosols and/or fogs during winter inversions is a likely reason for the relatively large PG in winter. A seasonal dependence of the PG also makes it difficult to analyze a potential influence of the B_y polarity on the PG using anything other than spectral methods. Figure 8 shows the PSD of the normalized PG that is displayed in Figure 7. Figure 8 demonstrates that no significant peak around the period of solar equatorial rotation (∼27 days) was observed in the PG data. Therefore, unlike the studies in polar regions (Burns et al., 2008; Lam and Tinsley, 2016), no influence of the B_y polarity on the PG was observed at the considered middle-latitude stations. However, a local, nondominant peak at the period of around ∼25.7 days corresponding to solar rotation at high solar latitudes was found in the PSD for all three stations. It should be remembered that a peak at this period was also observed in the PSD of v_x. However, the spectra of v_x and the PG do not exhibit similarities at other periods.

FIGURE 7

FIGURE 7. Normalized PG (proxy for fair weather) at Studenec (A), Panska Ves (B), and LS (C).

FIGURE 8

FIGURE 8. Power spectral densities (PSD) of the normalized PG (proxy for fair weather) measured at Studenec (A), Panská Ves (B) and on LS (C).

The PSD of PG data measured in Panská Ves and Studenec also exhibits sharp peaks at the periods of ∼0.5 and 1 day. It means that the Carnegie curve has two peaks (Nicoll et al., 2019). Interestingly, the PSD of PG data recorded in Panská Ves (partly also in Studenec) has a peak at the period of 7 days. It could possibly indicate a weekend reduction of industrial aerosol concentrations (Silva et al., 2014) as Panská Ves lies the closest to urban areas of all these three stations. A detailed analysis of this phenomenon is outside the scope of the present study.

Phase Relations and Lightning Activity

As we have found a relatively good correspondence between the power spectral density peaks at the period of ∼27 days for B_x and B_y and for the SCR counts corresponding to the low-energy charged particles and neutrons measured using the SEVAN [1,0,0] and [0,1,0] combinations, respectively, and for the NM counts (compare Figure 3A with Figures 5A, B, D), it is useful to investigate phase relations (time delays) between these variables around the period of ∼27 days. Figure 9 presents the CSD of the normalized B_x and NM counts, calculated according to Eq 2. Only data from the time interval from June 10, 2016 to April 15, 2019 were used to analyze the same period of time for both data series. Figures 9A, B show the absolute value and the phase of the CSD, respectively. The phase difference between the NM and B_x counts fluctuates around 40° at the period of ∼27 days, which clearly dominates in the cross-spectrum. In other words, the count variations observed using the NM advance the variations of the near-Earth B_x component of the HMF at the period of ∼27 days by approximately 3 days. Very similar results were obtained for cross-spectral analyses of the SEVAN [1,0,0] and [0,1,0] combinations and B_x (not shown). Calculating the cross-spectrum of B_x and B_y, it is easy to verify that the phase difference between B_x and B_y is very close to 180°. Consequently, the count variations recorded using the NM and the SEVAN [1,0,0] and [0,1,0] combinations lag the variations of B_y. The phase difference is around –140° and exhibits a frequency dependence, as is demonstrated in Figures 9C–H, which show the absolute values and phase differences of the CSD for the NM–B_y, SEVAN [1,0,0]–B_y, and SEVAN [0,1,0]–B_y data series. Asterisks in the individual plots of Figure 9 mark the maximum of the absolute values of the CSD and values larger than ½ of this maximum, including the corresponding phase differences. The phase differences change around the maxima; nevertheless, they are relatively stable for the frequencies (periods) of interest, marked by the asterisks, unlike the values of phase differences outside the maxima.

FIGURE 9

FIGURE 9. Cross-spectral densities (absolute values and phase differences) of (A, B) normalized pressure-corrected NM counts and Bx, (C, D) normalized pressure-corrected NM counts and B_y, (E, F) normalized pressure-corrected SEVAN [1,0,0] counts and B_y, and (G, H) normalized pressure-corrected SEVAN [0,1,0] counts and B_y. Asterisks mark the peak of absolute values of CSD, including its flanks, and the corresponding phase difference.

It is useful to remember that it is assumed that SCRs might play an important role in lightning initiation by providing seed energetic electrons and the ionization necessary for breakdown processes (e.g., Dwyer and Uman, 2014). Recently, Shao et al. (2020) performed a detailed radio frequency interferometry of the fast positive breakdown process that started a normal intracloud lightning and concluded that the breakdown process was ignited by a cosmic ray shower. Thus, there is a possibility that the discovered ∼27-day modulation of SCRs that is primarily caused by the ∼27-day periodicity in the HMF and solar wind, which influences the primary CRs entering the Earth’s atmosphere, could also modulate lightning occurrence. Indeed, some previous statistical studies indicate that the polarity of B_y and corotating interaction regions that often embed the HCS modulate lightning activity in middle latitudes (Owens et al., 2014; Owens et al., 2015; Miyahara et al., 2018). It should be noted that the polarity of B_y (B_x) reverses in the HCS.

On the other hand, the high electric field in thunderclouds might also be responsible for enhancements of SCRs. Such enhancements were also observed on Lomnický Štít (Kudela et al., 2017; Chum et al., 2020). In these cases, however, the seed energetic electrons generated by cosmic rays were likely multiplied in the process called relativistic runaway electron avalanche (RREA) suggested by Gurevich et al. (1992). Thunderstorm ground enhancements (TGEs) of SCRs were usually observed only in SEVAN channel 1, except one extreme event, for which the TGE was also reliably detected using the NM and the SEVAN [0,1,0] combination, which indicated a possibility of photonuclear reactions in the atmosphere and/or the material surrounding the detectors (Chum et al., 2020). It should be remembered in this respect that the significant enhancements of SCRs due to high electric fields—the several minutes–long TGEs—represented outliers in our dataset and were removed from the analysis, as was described in the section “Methods.” Therefore, the observed ∼27-day periodicity of SCRs is modulated by the HMF/SW and not by thunderstorms.

We investigated the periodicity of lightning activity in the area with the coordinates 48–51°N, 11–23°E, covering roughly Czechia and Slovakia. Figure 10A presents the number of lightnings recorded by the WWLLN in subsequent 1-h intervals for the period 2016–2019. Obviously, the lightning counts are high in the summer seasons. The seasonal variability produces a spectral peak around the period of one year and partly around lower harmonics (∼half, ∼third, ∼quarter, and ∼fifth of the year). The spectral peaks related to the seasonal variability are therefore sufficiently far from the period of ∼27 days on the frequency (period) axis. To calculate the spectra, lightning counts were first normalized using the method described in the section “Methods.” Figure 10B shows the PSD of normalized lightning counts. The most distinct and sharp spectral peak is observed for the period of 1 day (Figure 10B) as the lightning activity usually peaked between 14–15 UT over Czechia and Slovakia in the analyzed period. Besides the dominant peak at the 1-day period, several relatively flat peaks of similar amplitudes also occurred. One of these flat peaks is at the period of ∼27 days. Figures 10C, D display results of cross-spectral analysis between the lightning counts and B_y; specifically, the absolute values of CSD and phase differences are shown, respectively. The distinct peak at the period of ∼27 days of absolute values of the CSD is mainly caused by the high power spectral density of B_y at this period. Asterisks in the individual plots of Figure 9 mark the maximum of the absolute values of the CSD at the period of ∼27 days and values larger than ½ of this maximum, including the corresponding phase differences. It should be noted that the phase difference around the maximum is relatively stable. The observed phase difference is about 180°, which indicates that a probability of lightning occurrence is in antiphase with B_y. In other words, the probability of lightning occurrence decreases for B_y > 0 and increases for B_y < 0. Consequently, it is in phase with B_x (not shown) as fluctuations of B_x and B_y are anticorrelated at the period of ∼27 days. It should be emphasized that the dependence of lightning activity on the values/polarity of B_x and B_y is only weak as the spectral peak at ∼27 days in the spectrum of lightning counts (Figure 10B) is not dominant. Nevertheless, this finding is an interesting result which is in agreement with the previous statistical study performed by Owens et al. (2014), who showed that the lightning occurrence over the UK was modulated by the polarity of the HMF and found that the lightning occurrence was lower for the orientation of the HMF away from the Sun (B_x < 0 and B_y > 0) than it was for the orientation toward the Sun (B_x > 0 and B_y < 0).

FIGURE 10

FIGURE 10. Analysis of lightning activity. (A) Lightning counts over Czechia and Slovakia. (B) Power spectral density (PSD) of normalized lightning counts. (C, D) Cross-spectral densities (absolute values and phase differences) of normalized lightning counts and By. Asterisks mark the peak of absolute values of CSD, including its flanks, and the corresponding phase difference.

Conclusion

Effects of the HMF on SCRs and on atmospheric electricity at middle latitudes were studied by using measurements of SCRs in Slovakia (LS), the electric field in Slovakia (LS) and Czechia (Panská Ves and Studenec), and lightning activity in Czechia and Slovakia. By analyzing the spectra of the near-Earth HMF, velocities of SW, atmospheric electric field, and SCRs measured using the NM and SEVAN detectors located on the mount of Lomnický Štít, Slovakia, it was found that the variations with the period of solar rotation (∼27 days) and its second harmonic (∼13.5 days) are well expressed both in the HMF and SW and in SCR data. Especially, count variations registered by the NM and by the SEVAN low-energy and neutron channels ([1,0,0] and [0,1,0] combinations, respectively) exhibit distinct peaks at the period of ∼27 days. We also investigated the phase difference between the B_x and B_y components of the HMF and the variations of the NM and SEVAN counts and found that the variations of NM and SEVAN [1,0,0] and [0,1,0] data advance, on average, the fluctuations of the B_x component of the HMF by roughly 40°, which corresponds to about a 3-day time-shift at the period of 27 days. As B_x and B_y are practically in the antiphase as predicted by Parker (1958), the variations of NM, SEVAN [1,0,0], and SEVAN [0,1,0] data lag the fluctuations of the B_y component of the HMF by roughly –140°.

No reliable signatures of the ∼27-day periodicity were found in the electric field (PG) measurements in Czechia and Slovakia. Therefore, we conclude that the B_y component of the HMF does not influence the PG measured at these middle-latitude stations. This differs from the observation at polar stations, mainly stations in the winter Antarctic, where the fair-weather PG substantially depended on the B_y polarity (Mansurov et al., 1974; Burns et al., 2008). On the other hand, minor signatures of ∼25.7-day periodicities were found in the analyzed PG observations. Such a periodicity was also observed in the SW velocity. Because of the large differences between the PG and SW spectra, it is difficult to reliably determine whether the 25.7-day periodicity in PG data is induced by SW, as it is likely a random coincidence. Our observations show that variations of the PG measured in Czechia and Slovakia are likely fully controlled by processes in the troposphere.

The lightning activity in Slovakia and Czechia and its periodicity were investigated using data from the WWLLN. We observed a minor peak of the lightning activity at the period of solar rotation of about 27 days. The fluctuation of lightning counts at this period was in phase with the B_x component of the HMF and was in antiphase with the B_y component. In other words, lightning was more probable for the orientation of the HMF toward the Sun (B_x > 0 and B_y < 0) than for the orientation away from the Sun (B_x < 0 and B_y > 0). It should be noted that this effect is weak and that the lightning activity is mainly controlled by processes in the troposphere. A similar relation was also observed for the lightning occurrence in the UK by Owens et al. (2014), who used, however, a statistical approach in their study. The fact that similar results were obtained using different methods and for different regions enhances the credibility of these findings.

A more detailed investigation of the phase differences between the measured SCR, lightning occurrence, and B_y (B_x) for different phases of the solar cycle, together with investigation of other characteristics of solar activity, and analysis of the physical mechanisms leading to the observed relations is a potential subject for future study.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: the SW and HMF data can be found at NASA/GSFC’s Space Physics Data Facility’s OMNIWeb (https://omniweb.gsfc.nasa.gov/). The SEVAN and NM data and PG data from Lomnický Štít can be downloaded from http://data.space.saske.sk/status/ (access is provided by R. Langer, langer@saske.sk, on request). PG data from Czechia are available at https://www.ufa.cas.cz/en/institute-structure/department-of-ionosphere-and-aeronomy/measurements-operated-by-department-of-ionosphere-and-aeronomy/.

Author Contributions

JC designed and wrote the manuscript and performed most of the analysis. MK, RL, and IS are responsible for and provided the SCR data. IK provided the lightning counts from the WWLLN. DS helped with solar wind data. JR helped with the wavelet transform. All authors read and approved the submitted version.

Funding

Support under the grants SAV-18-04 and SAV-21-01 by the Czech Academy of Sciences is acknowledged. IK also acknowledges the support under the project 20-09671S of the Czech Science Foundation.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are grateful to Jiří Baše for maintaining the electric field measurements. JC also acknowledges valuable discussions with Keri Nicoll that preceded this study.

Abbreviations

CR, cosmic rays; CSD, cross-spectral density; EFM, electric field mill; GSE, geocentric solar ecliptic (coordinate system); HCS, heliospheric current sheet; HMF, heliospheric magnetic field; LS, Lomnický Štít; NM, neutron monitor; PG, potential gradient; PSD, power spectral density; SCR, secondary cosmic rays; SW, solar wind; WWLLN, World Wide Lightning Location Network; WT, wavelet transform.

References

Bazilevskaya, G., Broomhall, A.-M., Elsworth, Y., and Nakariakov, V. M. (2014). A Combined Analysis of the Observational Aspects of the Quasi-Biennial Oscillation in Solar Magnetic Activity. Space Sci. Rev. 186, 359–386. doi:10.1007/s11214-014-0068-0

CrossRef Full Text | Google Scholar

Burns, G. B., Tinsley, B. A., French, W. J. R., Troshichev, O. A., and Frank-Kamenetsky, A. V. (2008). Atmospheric Circuit Influences on Ground-Level Pressure in the Antarctic and Arctic. J. Geophys. Res. 113, D15112. doi:10.1029/2007JD009618

CrossRef Full Text | Google Scholar

Chilingarian, A., Babayan, V., Karapetyan, T., Mailyan, B., Sargsyan, B., and Zazyan, M. (2018). The SEVAN Worldwide Network of Particle Detectors: 10 Years of Operation. Adv. Space Res. 61, 2680–2696. doi:10.1016/j.asr.2018.02.030

CrossRef Full Text | Google Scholar

Chilingarian, A., Hovsepyan, G., Karapetyan, G., and Zazyan, M. (2021). Stopping Muon Effect and Estimation of Intracloud Electric Field. Astroparticle Phys. 124, 102505. doi:10.1016/j.astropartphys.2020.102505

CrossRef Full Text | Google Scholar

Chilingarian, A., and Reymers, A. (2008). Investigations of the Response of Hybrid Particle Detectors for the Space Environmental Viewing and Analysis Network (SEVAN). Ann. Geophys. 26, 249–257. Available at: www.ann-geophys.net/26/249/2008/. doi:10.5194/angeo-26-249-2008

CrossRef Full Text | Google Scholar

Chowdhury, P., and Kudela, K. (2018). Quasi-periodicities in Cosmic Rays and Time Lag with the Solar Activity at a Middle Latitude Neutron Monitor: 1982-2017. Astrophys. Space Sci. 363, 250. doi:10.1007/s10509-018-3467-y

CrossRef Full Text | Google Scholar

Chum, J., Langer, R., Baše, J., Kollárik, M., Strhárský, I., Diendorfer, G., et al. (2020). Significant Enhancements of Secondary Cosmic Rays and Electric Field at the High Mountain Peak of Lomnický Štít in High Tatras during Thunderstorms. Earth Planets Space 72, 28. doi:10.1186/s40623-020-01155-9

CrossRef Full Text | Google Scholar

Dwyer, J. R., and Uman, M. A. (2014). The Physics of Lightning. Phys. Rep. 534 (4), 147–241. doi:10.1016/j.physrep.2013.09.004

CrossRef Full Text | Google Scholar

Gosling, J. T., and Pizzo, V. J. (1999). Formation and Evolution of Corotating Interaction Regions and Their Three Dimensional Structure. Space Sci. Rev. 89, 21–52. doi:10.1023/A:100529171190010.1007/978-94-017-1179-1_3

CrossRef Full Text | Google Scholar

Gosse, J., and Klein, J. (2015). “Terrestrial Cosmogenic Nuclide Dating,” in Encyclopedia of Scientific Dating Methods. Editors W. J. Rink, and J. W. Thompson (New York: Springer), 29, 799–813. doi:10.1007/978-94-007-6304-3_148

CrossRef Full Text | Google Scholar

Gray, L. J., Beer, J., Geller, M., Haigh, D. J., Lockwood, M., Matthes, K., et al. (2010). Solar Influences on Climate. Rev. Geophys. 48, RG4001. doi:10.1029/2009RG000282

CrossRef Full Text | Google Scholar

Grieder, P. (2001). Cosmic Rays at Earth, 1. Amsterdam: Elsevier.

Gurevich, A. V., Milikh, G. M., and Roussel-Dupre, R. A. (1992). Runaway Electron Mechanism of Air Breakdown and Preconditioning during a Thunderstorm. Phys. Lett. A 165 (5-6), 463–468. doi:10.1016/0375-9601(92)90348-P

CrossRef Full Text | Google Scholar

Kirkby, J. (2008). Cosmic Rays and Climate. Surv. Geophys. 28, 333–375. doi:10.1007/s10712-008-9030-610.1007/s10712-008-9030-6

CrossRef Full Text | Google Scholar

Kristjánsson, J. E., Stjern, C. W., Stordal, F., Fjæraa, A. M., Myhre, G., and Jónasson, K. (2008). Cosmic Rays, Cloud Condensation Nuclei and Clouds - a Reassessment Using MODIS Data. Atmos. Chem. Phys. 8 (24), 7373–7387. doi:10.5194/acp-8-7373-2008

CrossRef Full Text | Google Scholar

Kudela, K., Chum, J., Kollárik, M., Langer, R., Strhárský, I., and Baše, J. (2017). Correlations between Secondary Cosmic Ray Rates and Strong Electric Fields at Lomnický Štít. J. Geophys. Res. Atmos. 122, 700–710. doi:10.1002/2016JD026439

CrossRef Full Text | Google Scholar

Kudela, K., and Langer, R. (2009). Cosmic ray Measurements in High Tatra Mountains: 1957-2007. Adv. Space Res. 44 (10), 1166–1172. doi:10.1016/j.asr.2008.11.028

CrossRef Full Text | Google Scholar

Kudela, K., Mavromichalaki, H., Papaioannou, A., and Gerontidou, M. (2010). On Mid-term Periodicities in Cosmic Rays. Sol. Phys. 266, 173–180. doi:10.1007/s11207-010-9598-0

CrossRef Full Text | Google Scholar

Kudela, K., Rybák, J., Antalová, A., and Storini, M. (2002). Time Evolution of Low-Frequency Periodicities in Cosmic ray Intensity. Solar Phys. 205, 165–175. doi:10.1023/A:1013869322693

CrossRef Full Text | Google Scholar

Lam, M. M., and Tinsley, B. A. (2016). Solar Wind-Atmospheric Electricity-Cloud Microphysics Connections to Weather and Climate. J. Atmos. Solar-Terrestrial Phys. 149, 277–290. doi:10.1016/j.jastp.2015.10.019

CrossRef Full Text | Google Scholar

López-Comazzi, A., and Blanco, J. J. (2020). Short-Term Periodicities Observed in Neutron Monitor Counting Rates. Sol. Phys. 295, 81. doi:10.1007/s11207-020-01649-5

CrossRef Full Text | Google Scholar

Mansurov, S. M., Mansurova, L. G., Mansurov, G. S., Mikhnevich, V. V., and Visotsky, A. M. (1974). North–south Asymmetryofgeomagneticandtroposphericevents. J.Atmos.Terr.Phys. 36 (1), 1957–1962. doi:10.1016/0021-9169(74)90182-2

CrossRef Full Text | Google Scholar

Miyahara, H., Kataoka, R., Mikami, T., Zaiki, M., Hirano, J., Yoshimura, M., et al. (2018). Solar Rotational Cycle in Lightning Activity in Japan during the 18-19th Centuries. Ann. Geophys. 36, 633–640. doi:10.5194/angeo-36-633-2018

CrossRef Full Text | Google Scholar

Nicoll, K. A., Harrison, R. G., Barta, V., Bor, J., Brugge, R., Chillingarian, A., et al. (2019). A Global Atmospheric Electricity Monitoring Network for Climate and Geophysical Research. J. Atmos. Solar-Terrestrial Phys. 184, 18–29. doi:10.1016/j.jastp.2019.01.003

CrossRef Full Text | Google Scholar

Nicoll, K. A., and Harrison, R. G. (2009). Vertical Current Flow through Extensive Layer Clouds. J. Atmos. Solar-Terrestrial Phys. 71 (12), 2040–2046. doi:10.1016/j.jastp.2009.09.011

CrossRef Full Text | Google Scholar

Owens, M. J., and Forsyth, R. J. (2013). The Heliospheric Magnetic Field. Living Rev. Solar Phys. 10, 5. doi:10.12942/lrsp-2013-5

CrossRef Full Text | Google Scholar

Owens, M. J., Scott, C. J., Bennett, A. J., Thomas, S. R., Lockwood, M., Harrison, R. G., et al. (2015). Lightning as a Space-Weather hazard: UK Thunderstorm Activity Modulated by the Passage of the Heliospheric Current Sheet. Geophys. Res. Lett. 42, 9624–9632. doi:10.1002/2015GL066802

CrossRef Full Text | Google Scholar

Owens, M. J., Scott, C. J., Lockwood, M., Barnard, L., Harrison, R. G., Nicoll, K., et al. (2014). Modulation of UK Lightning by Heliospheric Magnetic Field Polarity. Environ. Res. Lett. 9 (11), 115009. doi:10.1088/1748-9326/9/11/115009

CrossRef Full Text | Google Scholar

Parker, E. N. (1958). Dynamics of the Interplanetary Gas and Magnetic Fields. ApJ 128, 664–676. doi:10.1086/146579

CrossRef Full Text | Google Scholar

Riádigos, I., García‐Castro, D., González‐Díaz, D., and Pérez‐Muñuzuri, V. (2020). Atmospheric Temperature Effect in Secondary Cosmic Rays Observed with a 2 M 2 Ground‐Based tRPC Detector. Earth Space Sci. 7, e2020EA001131. doi:10.1029/2020EA001131

CrossRef Full Text | Google Scholar

Rieger, E., Share, G. H., Forrest, D. J., Kanbach, G., Reppin, C., and Chupp, E. L. (1984). A 154-day Periodicity in the Occurrence of Hard Solar Flares? Nature 312, 623–625. doi:10.1038/312623a0

CrossRef Full Text | Google Scholar

Rodger, C. J., Brundell, J. B., Dowden, R. L., and Thomson, N. R. (2004). Location Accuracy of Long Distance VLF Lightning Locationnetwork. Ann. Geophys. 22 (3), 747–758. doi:10.5194/angeo-22-747-2004

CrossRef Full Text | Google Scholar

Rycroft, M. J., Israelsson, S., and Price, C. (2000). The Global Atmospheric Electric Circuit, Solar Activity and Climate Change. J. Atmos. Solar-Terrestrial Phys. 62, 1563–1576. doi:10.1016/S1364-6826(00)00112-7

CrossRef Full Text | Google Scholar

Saad Farid, A. I. (2019). High Frequency Spectral Features of Galactic Cosmic Rays at Different Rigidities during the Ascending and Maximum Phases of the Solar Cycle 24. Astrophys Space Sci. 364, 57. doi:10.1007/s10509-019-3544-x

CrossRef Full Text | Google Scholar

Shao, X. M., Ho, C., Bowers, G., Blaine, W., and Dingus, B. (2020). Lightning Interferometry Uncertainty, Beam Steering Interferometry, and Evidence of Lightning Being Ignited by a Cosmic ray Shower. J. Geophys. Res. Atmos. 125, e2019JD032273. doi:10.1029/2019JD032273

CrossRef Full Text | Google Scholar

Shea, M. A., and Smart, D. F. (2000). Fifty Years of Cosmic Radiation Data. Space Sci. Rev. 93, 229–262. doi:10.1023/A:1026500713452

CrossRef Full Text | Google Scholar

Silva, H. G., Conceição, R., Melgão, M., Nicoll, K., Mendes, P. B., Tlemçani, M., et al. (2014). Atmospheric Electric Field Measurements in Urban Environment and the Pollutant Aerosol Weekly Dependence. Environ. Res. Lett. 9 (11), 114025. doi:10.1088/1748-9326/9/11/114025

CrossRef Full Text | Google Scholar

Svensmark, H., Bondo, T., and Svensmark, J. (2009). Cosmic ray Decreases Affect Atmospheric Aerosols and Clouds. Geophys. Res. Lett. 36, L15101. doi:10.1029/2009GL038429

CrossRef Full Text | Google Scholar

Usoskin, I. G., Kananen, H., Mursula, K., Tanskanen, P., and Kovaltsov, G. A. (1998). Correlative Study of Solar Activity and Cosmic ray Intensity. J. Geophys. Res. 103 (A5), 9567–9574. doi:10.1029/97JA03782

CrossRef Full Text | Google Scholar

Voiculescu, M., Usoskin, I., and Condurache-Bota, S. (2013). Clouds Blown by the Solar Wind. Environ. Res. Lett. 8, 045032. doi:10.1088/1748-9326/8/4/045032

CrossRef Full Text | Google Scholar

Voiculescu, M., and Usoskin, I. (2012). Persistent Solar Signatures in Cloud Cover: Spatial and Temporal Analysis. Environ. Res. Lett. 7, 044004. doi:10.1088/1748-9326/7/4/044004

CrossRef Full Text | Google Scholar

Yang, Z.-Y., and Sheu, R.-J. (2020). An In-Depth Analysis of Aviation Route Doses for the Longest Distance Flight from Taiwan. Radiat. Phys. Chem. 168, 108548. doi:10.1016/j.radphyschem.2019.108548

CrossRef Full Text | Google Scholar