On the Geomagnetic Field Line Resonance Eigenfrequency Variations during Seismic Event

: In this paper, we report high statistical evidence for a seismo–ionosphere effects occurring in conjunction with an earthquake. This ﬁnding supports a lithosphere-magnetosphere coupling mechanism producing a plasma density variation along the magnetic ﬁeld lines, mechanically produced by atmospheric acoustic gravity waves (AGWs) impinging the ionosphere. We have analysed a large sample of earthquakes (EQ) using ground magnetometers data: in 28 of 42 analysed case events, we detect a temporary stepwise decrease ( ∆ f ) of the magnetospheric ﬁeld line resonance (FLR) eigenfrequency ( f ∗ ). ∆ f decreases of ∼ 5–25 mHz during ∼ 20–35 min following the time of the EQ. We present an analytical model for f ∗ , able to reproduce the behaviour observed during the EQ. Our work is in agreement with recent results conﬁrming co-seismic direct coupling between lithosphere, ionosphere and magnetosphere opening the way to new remote sensing methods, from space/ground, of the earth seismic activity.


Introduction
The study of the physical process connected to the preparation and onset of an earthquake is a topic of increasing interest among the scientific community, also in view of the societal impact of these phenomena. One of the challenges of these studies is to identify physical phenomena which can be directly connected, without ambiguity, with the earthquake geographical location and time window. Most of the evidence in the literature is, indeed, of a statistical nature, while event based, causal observations of the connection among ground, ionosphere and magnetosphere are much more difficult to be convincingly demonstrated. Regarding the statistical evidence, one of the the most interesting and promising result is related to electromagnetic and ionospheric disturbances occurring before and during seismic activities. Examples of these results are the experimental investigation of the lithosphere-ionosphere-magnetosphere coupling [1][2][3] with the observation of "anomalous" pulses of electromagnetic (EM) emissions in the frequency interval between a few Hz and up to few tens of kHz, as well as the more recent observations of changes of the density of the charged trapped particles registered by satellites [4,5]. More recently, investigations of earthquake preparation phenomena using data registered by the DEMETER satellite provided statistical evidence for spectral damping of VLF (very low frequency) radio signals at F-region altitudes and within a radius of 1000-5000 km from the earthquake epicenter, about 0-3 weeks before the event [6]. In addition, using DEMETER electric and magnetic field data, Bertello et al. [7] found an EM wave at ∼300 Hz propagating 2 days before the L'Aquila 2009 earthquake event. In order to better understand the physical processes present in the lithosphere-atmosphere-magnetosphere interactions, many studies focused on the disturbances induced on the atmospheric electric field [8], on the anomalous geomagnetic pulsations [9,10], as well as on other anomalous disturbances in the ionosphere [11] and magnetosphere [12,13]. Searches for possible seismo-ionospheric effects were also performed before the earthquake by the satellite Interkosmos-19 operating at F-region altitudes, providing for the first time evidence for an increase of both the intensity of VLF noise, in the frequency range between 140 Hz to 15 kHz, and for disturbances of the electron density at a distance from the epicentres up to a few 1000 km [3,14]. More recently, Carbone et al. [15], Piersanti et al. [16] started the development of an analytical model of the coupling between lithosphere-atmosphere-ionosphere-magnetosphere to be submitted to detailed experimental verification (M.I.L.C.). In the M.I.L.C. model the coupling during active seismic conditions is described by the onset of atmospheric and ionospheric EM and particle anomalies: a first successful test of the model was the analysis of the 2018 Bayan EQ, when a series of correlated phenomena were detected both by ground sensors and by low earth orbiting satellites (∼500 km) around the time of EQ occurrence. The authors explained and modelled the experimental observations as due to the generation of an acoustic gravity wave (AGW) induced by the EQ which mechanically perturbed the ionospheric medium causing both an EM emissions and plasma waves. Interestingly, the model predicts a clear decrease of the magnetospheric FLR f * in concomitance of the EQ occurrence, which has also been observed. This phenomenon was never reported before in the literature and it is particularly interesting, since it represents a direct, unambiguous evidence of the connection between the lithosphere and the magnetosphere, which can be used both for the analysis of coseismic as well as of precursor phenomena. Following the result on the 2018 Bayan EQ, we started a systematic study of this phenomena using 42 EQ in the time span from 2001-07-17 and 2020-08-31. This paper presents the result of this study, in which we analyze the f * variations using ground magnetometers observations, and we explain the results of these experimental observations with an analytical model describing the f * behaviour during active seismic conditions.

Data and Methods
The ground magnetometers used for the present analysis come from both INTER-MAGNET and SUPERMAG magnetometer array networks, which are consortium of observatories guaranteeing a common standard data release to the scientific community, leading to possible comparison among measurements at different observation points. In our analysis we have used 1 s time resolution data.
To evaluate the f * , we studied the cross-space spectrum [17] between the North-South magnetic field components observed at two geomagnetic observatories close enough to the EQ epicenter location (see Table 1). It is well known that, at the eigenfrequency of a field line centered between two neighboring stations having almost the same magnetic longitude, the phase difference maximizes [17,18]. Waters et al. [19] showed that the patterns of the maximum phase differences in the cross-phase spectrograms were observed consistently from day to day in the dayside region over baselines of about 100 km in the magnetic meridian. Green et al. [18] also reported that, among the several methods that determine the resonant frequency, the phase shift is least affected by geologic inhomogeneity and consistently defines the resonant frequencies.

FLR Frequency Behaviour during Seismic Events
We have evaluated the FLR frequency behaviour for 42 low latitudes EQs (below 39 • of geographical latitude) in the time span from 2001-07-17 and 2020-08-31. A part of these EQ (first 32 events) belongs to the sample selected by Battiston et al. [4] using POES satellite data. All the EQs have been chosen as result of a cross-check with the planetary geomagnetic K p index [20] in order to exclude any possible f * variation of solar origin. Table 1 summarizes the results. First of all, the K p index ranges between 0 and 2, indicating that any possible f * variation of solar origin can be reasonably neglected. The EQ are characterized by a magnitude (M) greater than 5. Then, we found 28 cases out of 42 in which there is a clear variation of the estimated f * (indicated with the X). In 4 cases it was not possible to correctly evaluate f * because of the post-sunset occurrence of the EQ (indicated with the N A). In fact, as explained in Menk et al. [21], the determination of the FLR frequency usually fails in the nightside regions. Finally, no f * variations has been detected for 10 case events (indicated with −). Figures 1 and 2 show four examples of FLR eigenfrequency time dependence in a time window around the EQ occurrence (red dashed line) using the cross-phase spectrogram. Colours are representative of the phase difference between the two stations selected for the f * evaluation.  Each spectrum has been evaluated over one hour interval. The spectra have been smoothed both in time and frequency domains (7 frequency bands and 15 temporal bands). The red vertical line represents the earthquake occurrence time. The green dashed line shows the occurrence of f * decrease ∼2 h before the EQ main shock. The top caption reports the INTERMAGNET ground station codices used for the evaluation of the dynamical cross-phase spectrogram. The color-bar represents the phase difference in degrees between the equatorward and poleward ground magnetometer.
As expected, at low latitudes the FLR eigenfrequency is around 110 mHz [19,22,23]. For each event, 2 ± 1 min after the earthquake occurrence (red dashed vertical line) there is a clear decrease of f * . In fact, the upper (a), the middle (b) and the lower panels (c) show variations of ∼−10 mHz, ∼−25 mHz and ∼−12 mHz, respectively. The time duration of such variation is ∼15 min for the first (panel a) and the third (panel c) event, and ∼30 min for the second event (panel b). It is worth highlighting here that in 97% the FLR frequency variation were characterized by a single decrease coincident with the EQ occurrence, while in the remaining 3% was featured by a double f * reduction as reported in Figure 2. In fact, in addition to the decrease of the FLR eigenfrequency at the moment of the EQ occurrence (red dashed vertical line), a clear reduction of f * is also visible less than two hours before (green vertical dashed line). The variation is of ∼−10 mHz both for the coseismic f * decrease as well as for the precursor f * decrease, while the time duration is ∼40 min for the precursor phenomenon and ∼25 min for the coseismic phenomenon. Figure 3 shows the statistical analysis of the EQ events characterized by a FLR decrease in terms of frequency variation (δ f ) and relative time duration (δT). It can be easily seen that the typical δT of the eigenfrequency decrease (panel a) is between 25 and 35 min. On the other hand, δ f on average shows variations of ∼10 mHz. Figure 3c) shows the probability density (dP) of δ f as a function of δT. dP has been estimated constructing bivariate histograms and using a kernel density estimator (e.g., [24]) with the following bin sizes: δ f , 3 mHz; δT, 3 min. It results that the co-seismic FLR eigenfrequency variation is characterized by a frequency decrease of 12 ± 3 mHz and a time duration of 36 ± 3 min.

Discussion
In order to provide a quantitative explanation to the observed FLR eigenfrequency variations, we modelled the f * behaviour during the occurrence of an EQ. It is well known that a geomagnetic field line, with both ends fixed in the ionosphere can be sketched as a string whose frequency depends on both the magnetic field geometry and the plasma density along the field line [19,21,25,26]. Following the approach of Singer et al. [27], we have evaluated f * for an arbitrary magnetic field geometry, starting from the Magnetohydrodynamic (MHD) equations related to a stationary EM wave. By referring to the model reported in Appendix A, we have numerically solved the Equation (A8) using the magnetic equator as reference point V A (s 0 ) = V A (eq) (V A being the Alfvén speed). We have used the IGRF (International Geomagnetic Reference Field) model [28] for the internal Earth's magnetic field, the T01 model [29,30] for the external part of the Earth's magnetic field and a radial power law dependence for the plasma mass density, ρ/ρ eq = (r/r eq ) −3 [26]. The boundary conditions of fixed footpoints have been established at some level in the ionosphere, i.e., at h = 120 km altitude corresponding to the E-layer, where the Alfvén wave is assumed to be perfectly reflected [17]. Finally, the values of the eigenfrequency f have been obtained through Equation (A9) of Appendix A. Figure 4 shows the modelled diurnal eigenfrequency behaviour of a field line footprinted at λ mag = 20 • (λ mag being the magnetic latitude). Around noon, we modified the plasma density at the footprint of the field line using a gradient pressure (∇p den ) which produces a density variation of 15% lasting for about 10 min. Such ∇p den is the result of the application of the M.I.L.C. model to an EQ characterized by a magnitude M EQ = 6.5, a PGA = 0.6 g and a ∆t = 20 s (PGA and ∆t being the Peak Ground Acceleration and time duration of the EQ, respectively). The M.I.L.C. model is based on the assumption that an EQ creates an acoustic gravity wave, which propagates through the atmosphere. The pressure gradient induced by the AGW causes local instability in the ionospheric plasma density distribution, giving rise to both plasma and EM waves propagating up to the magnetosphere. In general, it is well known that the concurring contribution of the EM wave energy and/or of the plasma density variation produces a change in the local FLR eigen-frequency [19,31,32]. Indeed, Figure 4 shows that our model is able to produce a clear f * collapse in correspondence to the plasma density gradient. This result is the consequence of Equation (A9) according to which any variation in the local magnetic field and/or in the local plasma density produces a corresponding change in f * . It is important to remind here that the use of the IGRF and the T01 model in solving Equation (A8) produces a 5% maximum error in the evaluation of f * [23]. The result displayed in Figure 4 is consistent with the experimental FLR eigenfrequency behaviour detected in correspondence of an EQ. In fact, both Figure 1 and 2 show a sudden decrease of f * of ∼10 mHz, 2 ± 1 min after the EQ occurrence lasting for 20-30 min. Such result completely agree with the probability density bi-variate distribution in Figure 3c). However, we need to stress here that at low magnetic latitudes (0 • ≤ λ mag ≤ 30 • ) the geomagnetic field line is almost completely surrounded by the ionosphere. As a consequence any alteration in the ionospheric plasma density induces a variation in the corresponding eigenfrequency. Consequently, we do interpret the f * changes observed in our 28 EQ events (see Table 1) as caused by the ionospheric plasma density variation induced by the emission of a co-seismic AGW leading to a pressure gradient [15].
Finally, in the case of the absence of a co-seismic AGW emission, no possible f * variations can be detected (10 case event, see Table 1). Such a hypothesis is confirmed by Carbone et al. [15], showing that the atmospheric fluctuations excited by a generic seismic event on the top of the first layer of the atmosphere can be evanescent. In fact, depending on the characteristic parameters of the EQ (length of the fault, peak ground acceleration strong time duration and so on), a the propagation of the AGW up to the ionosphere can be prevented. In order to confirm such hypothesis, for these events, we analyzed the vertical atmospheric temperature profiles using the approach described in Piersanti et al. [16] to catch for possible AGW injection. Here, we display the analysis of the 19/12/2006 Sumatra EQ, since the remaining nine case events show similar results. Figure 5a shows the atmospheric vertical temperature profile (T) as obtained from ERA5, which is the 5th generation atmospheric data set produced by the European Centre for Medium-Range Weather Forecasts [33]. The temperature fluctuations (T ), evaluated as the difference between T and its 2 km moving average, show the expected minimum and maximum at the tropopause (∼18 km) and the stratopause [34], respectively. A similar behaviour can be found in both the Brunt − Väisälä frequency (N 2 ) and the potential energy density (E P ) value [35]. The lack of any possible wave behaviour in T confirms the absence of AGW (and reference therein [36]) injected at the moment of the EQ occurrence. As a consequence, we can reasonably affirm that the missing of AGW prevents any possible variation of ionospheric plasma density distribution leading to the FLR eigenfrequency variation. Finally, it is worth noticing that the variation of f * does not show any dependence on earthquake magnitude (not shown). Such result agrees with Carbone et al. [15], who demonstrated that the emission of a non-evanescent AGW, generating FLR variation, does not depend on the individual earthquake parameter alone, but on both the combination of the length of the fault, the PGA, the time duration of the EQ, etc (see dispersion relation in Carbone et al. [15]), and the local atmospheric scale height.

Conclusions
In the last 20 years, many investigations focused on the possible identification of magnetospheric perturbations directly connected to earthquake occurrence ( [37] and reference therein). This paper presents the first evidence, via observation and modelling, of changes in magnetospheric FLR eigenfrequency associated to the EQ occurrence, demonstrating a causal connection between seismic phenomena and space based observables . We have analyzed more than 40 low latitudes EQ from 2000 to 2020, during quiet solar condition in order to search for magnetospheric signal associate to seismic activity. In 28 events, we found a clear sudden decrease of the magnetospheric FLR eigenfrequencies, while in 10 cases we did not find any f * variation. The proposed explanation is that the plasma density at the footprint of the field line magnetically connected to the EQ location was modified of about ∆ρ 15%, by a gradient pressure fluctuations ∇p den induced by a propagation of an AGW emitted during the EQ occurrence [15,16]. At low latitudes the magnetospheric field lines are fully surrounded by the ionosphere and the FLR eigenfrequencies, depending on both the magnetic field and the plasma density along the field line [19,31], is expected to decrease [23,26,32]. On the other hand, the possible explanation of the null f * variation observed in 10 case events, can be found in the lack of the vertical propagation of the AGW (evanescent) up to the ionosphere, as predicted by the Carbone et al. [15] analytical model.
In is interesting to note that the FLR decrease observed in one case some hours before the EQ occurrence ( Figure 2) could be due to various reasons, such as high-level seismic activity (especially for events characterized by a sequence of foreshocks before the main shocks), or to the outflow of radioactive gases (e.g., due to radon decay) by the Earth's surface [37,38]. Indeed, both these phenomena would be able to generate changes in atmospheric temperature and hence AGW formation (e.g., [39]). A similar result was found in Piersanti et al. [16] who found a decrease of f * 5 h before the EQ occurrence. They explained their observation in terms of the M.I.L.C. model, pointing out that any AGW can produce a variation of the ionospheric plasma density distribution (such as travelling ionospheric disturbances [40]) which in turns changes the Alfvén velocity along the field line giving rise to a change of the FLR eigenfrequency [19].
In conclusion, our results confirm analytically the direct coupling among lithosphere, ionosphere and magnetosphere during active seismic conditions, supporting the models introduced by [16,37] and opening the way to new remote sensing methods combining space and ground sensing of the earth seismic activity.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Field Line Resonance Eigenfrequency Equation
The external perturbation pressure δP AGW produced by the AGW, impinging on the ionospheric plasma from below, is able to modify the density field thus inducing a dynamics of plasma which, in turns, produces a perturbed magnetic field b = B − B 0 of the background magnetic field B 0 . As an order of magnitude estimate, from the Ampere law |∇b| ∼ β∇δρ, being δρ the fluctuating density and β is the plasma parameter. As a consequence, by using a relation between density and pressure fluctuation, for example the adiabatic gas law, it can be found |∇b| ∼ (βρ 2 0 /γP 0 B 0 )∇δP AGW , being ρ 0 and P 0 some reference values for density and pressure, and γ is the adiabatic index.
The corresponding low-frequency dynamics of the plasma can be roughly described by the ideal dissipationless MHD equations where v is the fluctuation velocity, j the current density, E the electric field fluctuations and P the internal pressure. The electric field fluctuations can be obtained by the Ohm's law by neglecting the Hall term and the electron pressure gradient because we are interested at the low-frequency evolution of plasma corresponding to scales much greater than the Larmor radius, so that E = −v × B. Furthermore, we are dealing with a low-β plasma [41], so that dynamical processes occurring within the ionospheric plasma cannot significant alter the background magnetic field, so that the internal pressure gradient ∇P can be neglected in the momentum Equation (A1). It is worthwhile to note that the same approximations are usually used to describe the plasma dynamics of low-β laboratory plasma, for example confined in Reversed Field Pinch devices (e.g., [42]).
To model the eigenfrequencies f * and amplitudes of low-frequency transverse waves, we use a linear model from Equation (A1), which describes the linear dynamics of fluctuations, namely where we used the Ampere's law ∇ × B = µ 0 j. Note that, by considering the plasma dynamics generated by AGW, the perturbed magnetic field can be viewed as generated by a small displacement ξ of the plasma [43], not by a compression of the field lines, so that b lye along the field line and produces the force in the momentum equation. The linear model (A2) neglects the background current density j 0 related to the background magnetic field [29,30]. In fact, as an order of magnitude estimate, the background current density |j 0 | 1.8 × 10 −10 A/m 2 results ten times lower than the current density |j| ∼ b/l 10 −9 A/m 2 (l is the scale length along the fluctuating magnetic field direction). Let us consider a geometry where the z-axis is directed between the field lines. However, both footpoints of a field line are fixed in the ionosphere, so that the displacement and the separation of adjacent lines can be described by a function h(x, y, z). From the last Equation (A2), using v = ∂ξ/∂t along the z-direction, we get the perturbed magnetic field (being e z the unit vector in the direction z) apart from a constant which can be cast to zero without losing generality. Then, from the momentum equation we obtain a relation for the displacement along the z direction By introducing the Alfvén velocity V A = B 0 / √ ρ 0 µ 0 and the normalized displacement ξ = ξ/h we finally obtain the wave equation for the displacement where e 0 represents the unit vector along the background magnetic field. If we consider now the ansatz where ξ behaves as e iωt , under the hypothesis that the field curvature is smooth enough so the function h is slowly variable, from Equation (A5) we obtain (e 0 · ∇) (e 0 · ∇)ξ + (e 0 · ∇) ln h 2 (e 0 · ∇)ξ + Introducing the coordinate s along the field line (e 0 · ∇) = −1 ∂/∂s, where is the characteristic length of the field line between two ionospheric footpoints, say using the coordinate system where 0 ≤ s ≤ 1, we obtain the characteristic value wave equation ∂ 2 ξ ∂s 2 + P(s) ∂ξ ∂s being P(s) = ∂ ln h 2 /∂s, a unknown function which depends on the coordinate along the field line. The characteristic frequencies f we are looking for, correspond to the characteristic values ω, which can be obtained once the Sturm-Liouville Equation (A7) is solved supplied by appropriate boundary conditions, for example the condition of fixed footpoints ξ (s ) = (∂ξ /∂s) s=s = 0 at both s = 0, 1.
To solve the characteristic value equation we can introduce a unknown eigenvalue λ by modifying the equation as where V A (s 0 ) is the value of the Alfvén speed in a point s 0 . The solution of Equation (A8) gives us the eigenvalue λ compatible with both the boundary conditions and a fixed value of V A (s 0 ). Finally, by a comparison of (A7) and (A8), the characteristic frequencies results to be It can be possible to analytically solve Equation (A8) for some particular geometries, by making explicit the function P(s) and estimate the value of V A (s 0 ). For example in a dipole field the azimuthal field line displacement, is proportional to h r sin θ [27], corresponding to a toroidal mode. However, we aimed to a direct comparison with real observations of the eigenfrequencies f , and this necessarily requires a numerical integration of Equation (A8), because we need the exact knowledge of the function P(s).