**3. 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 *Kp* index [20] in order to exclude any possible *f* ∗ variation of solar origin. Table 1 summarizes the results. First of all, the *Kp* 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 *NA*). 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.

**Figure 1.** The cross-phase dynamical spectrograms between two low-latitude ground stations near the earthquake epicenter: panel (**a**) Sumatra 16 September 2009 EQ; panel (**b**) Indonesia 6 January 2009 EQ; panel (**c**) Philippines 11 November 2002 EQ. Each spectrum has been evaluated over a 1 h interval. 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. In each panel 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.

**Figure 2.** The cross-phase dynamical spectrograms for the Indonesia 11 September 2008 earthquake. 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* <sup>∗</sup> decrease as well as for the precursor *f* <sup>∗</sup> 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 (*d*P) of *δ f* as a function of *δT*. *d*P 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.

**Figure 3.** Eigenfrequency variation (panel (**a**)) and relative time duration (panel (**b**)) data distributions of the FLR. Panel (**c**) shows the probability density distribution of the eigenfrequency variation as a function of the relative time duration.

#### **4. 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 *VA*(*s*0) = *VA*(*eq*) (*VA* 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*/*req*)−<sup>3</sup> [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 (∇*pden*) which produces a density variation of 15% lasting for about 10 min. Such ∇*pden* is the result of the application of the M.I.L.C. model to an EQ characterized by a magnitude *MEQ* = 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].

**Figure 4.** As an example, in the Figure, is shown the simulated behaviour of the FLR eigen-frequency obtained by our modeled for a potential EQ occurred at 20◦ magnetic latitude.

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* <sup>∗</sup> 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* − *Vais* ¨ *al*¨ *<sup>a</sup>*¨ frequency (*N*2) and the potential energy density (*EP*) 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.

**Figure 5.** Example of AGW analysis. Vertical profiles of: (**a**) temperature; (**b**) background temperature; (**c**) temperature deviation, (**d**) square term of *Brunt* − *Vais* ¨ *al*¨ *a*¨ frequency, and (**e**) potential energy at 13:00 UT on 19 December 2006.

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.

#### **5. 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 ∇*pden* 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.

**Author Contributions:** M.P. writing—original draft preparation, formal analysis, methodology, investigation and Conceptualization; W.J.B. writing—review and editing, methodology, investigation; V.C. writing—review and editing, validation; R.B. writing—review and editing, methodology and coordination; R.I. data curation and validation; P.U. writing—review and editing and project administration. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Ground geomagnetic observatory data are provided by INTER-MAGNET (www.intermagnet.org, accessed on 16 July 2021), and SUPERMAG (https://supermag. jhuapl.edu), accessed on 16 July 2021; Earthquake parameters comes from USGS data catalog (https://earthquake.usgs.gov/), accessed on 16 July 2021.

**Acknowledgments:** The authors thank the national institutes that support INTERMAGNET for promoting high standards of the magnetic observatory practice (www.intermagnet.org, accessed on 14 July 2021). The authors thank the provision of magnetometer data from SUPERMAG (https://supermag.jhuapl.edu, accessed on 14 July 2021). The parameters of the earthquake are provided by USGS data catalog (https://earthquake.usgs.gov/, accessed on 14 July 2021). M. Piersanti, R. Battiston and R. Iuppa thank the Italian Space Agency for the financial support under the contract ASI "LIMADOU Scienza+" n◦ 2020-31-HH.0. M. Piersanti thanks the ISSI-BJ project "The electromagnetic data validation and scientific application research based on CSES satellite" and Dragon 5 cooperation 2020-2024 (ID. 59236).

**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 *δPAGW* 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**<sup>0</sup> 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**| ∼ (*βρ*<sup>2</sup> <sup>0</sup>/*γP*0*B*0)∇*δPAGW*, being *ρ*<sup>0</sup> and *P*<sup>0</sup> 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

$$
\rho \left[ \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} \right]\_{\mathbf{v}} = \mathbf{j} \times \mathbf{B} - \nabla P
$$

$$
\frac{\partial \mathbf{B}}{\partial t}\_{\mathbf{j}} = \nabla \times \mathbf{E} \tag{A1}
$$

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

$$\begin{aligned} \rho\_0 \mu\_0 \frac{\partial \mathbf{v}}{\partial t} & \quad \simeq \quad (\nabla \times \mathbf{b}) \times \mathbf{B}\_0\\ \frac{\partial \mathbf{b}}{\partial t} & \quad \simeq \quad \nabla \times (\mathbf{v} \times \mathbf{B}\_0) \end{aligned} \tag{A2}$$

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**<sup>0</sup> related to the background magnetic field [29,30]. In fact, as an order of magnitude estimate, the background current density <sup>|</sup>**j**0<sup>|</sup> 1.8 <sup>×</sup> <sup>10</sup>−<sup>10</sup> A/m<sup>2</sup> results ten times lower than the current density <sup>|</sup>**j**| ∼ *<sup>b</sup>*/*l* <sup>10</sup>−<sup>9</sup> A/m2 (*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

$$\mathbf{b} = \nabla \times (\mathbf{j}\_{\varepsilon} \mathbf{e}\_{z} \times \mathbf{B}\_{0}) \tag{A3}$$

(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

$$
\rho\_0 \mu\_0 \frac{\partial^2 \xi}{\partial t^2} = (\mathbf{B}\_0 \cdot \nabla) [(\mathbf{B}\_0 \cdot \nabla)\xi] \tag{A4}
$$

By introducing the Alfvén velocity **<sup>V</sup>***<sup>A</sup>* <sup>=</sup> **<sup>B</sup>**0/√*ρ*0*μ*<sup>0</sup> and the normalized displacement *ξ* = *ξ*/*h* we finally obtain the wave equation for the displacement

$$\frac{\partial^2 \xi'}{\partial t^2} = V\_A^2(\mathbf{e}\_0 \cdot \nabla) \left[ h(\mathbf{e}\_0 \cdot \nabla) \xi' + \xi'(\mathbf{e}\_0 \cdot \nabla) h \right] \tag{A5}$$

where **e**<sup>0</sup> represents the unit vector along the background magnetic field.

If we consider now the ansatz where *ξ* behaves as *eiω<sup>t</sup>* , under the hypothesis that the field curvature is smooth enough so the function *h* is slowly variable, from Equation (A5) we obtain

$$\left[ (\mathbf{e}\_0 \cdot \nabla) \left[ (\mathbf{e}\_0 \cdot \nabla) \xi' \right] + \left[ (\mathbf{e}\_0 \cdot \nabla) \ln h^2 \right] \left[ (\mathbf{e}\_0 \cdot \nabla) \xi' \right] + \frac{\omega^2}{V\_A^2} \xi' \simeq 0 \tag{A6}$$

Introducing the coordinate *s* along the field line (**e**<sup>0</sup> · ∇) = -<sup>−</sup>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

$$\frac{\partial^2 \xi'}{\partial s^2} + P(s) \frac{\partial \xi'}{\partial s} + \left(\frac{\omega^2 \ell^2}{V\_A^2}\right) \xi' = 0 \tag{A7}$$

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

$$\frac{\partial^2 \xi'}{\partial s^2} + P(s) \frac{\partial \xi'}{\partial s} + \lambda \left[ \frac{V\_A^2(s\_0)}{V\_A^2(s)} \right] \xi' = 0 \tag{A8}$$

where *VA*(*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 *VA*(*s*0). Finally, by a comparison of (A7) and (A8), the characteristic frequencies results to be

$$f^\star = \frac{V\_A^2(s\_0)\sqrt{\lambda}}{\ell} \tag{A9}$$

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 *VA*(*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*).
