## 2 Materials and Methods
### Methods of Analysis  
Firstly, the _Swarm_ satellite data are pre-processed by MAgnetic _Swarm_ anomaly detection by Spline analysis (MASS) or electron density (_Ne_) Anomaly Detection (NeAD) algorithms to search for anomalies. The approach applied to the magnetic field and electron density data is basically the same, and only some parameters are settled differently to consider the different sample frequencies and signal-to-noise ratios properly. MASS and NeAD perform the following operations while taking into account the data along the track between \(-\)50\({}^{\circ}\) and \(+\)50\({}^{\circ}\) geomagnetic latitude (see Figure S1 in the Supplementary Materials as an example of data processing):  
1. First differences of the data: the difference between two consecutive data divided by the time interval between them (as a first order approximation of the temporal derivative). 2. Removal of the residual trend: a cubic spline coming from the fitting of the derivative of the data (as obtained from Step 1) is removed. 3. Analysis of the residuals obtained from Step 1 and Step 2 along the track: moving windows of 7\({}^{\circ}\) latitude are used, with an incremental shift of 1/5 of the window length (i.e., \(\sim\)1.4\({}^{\circ}\) of latitude). 4. The root mean square (rms) is calculated inside the window and compared with the whole track's Root Mean Square (RMS). The anomaly is defined if \(rms>k_{t}\cdot RMS\) (\(k_{t}\) is chosen as 2.5 or 3 if no-frequency investigation is performed; more details are given later on). In addition, a new frequency analysis is carried out, and the Fast Fourier Transform (FFT) is calculated for each moving window and for the whole spectrum. Three sub frequency bands are divided according to their period ranges: (2-10) s, (10-25) s, and (25-50) s. In each frequency band, the anomaly is defined if the mean intensity overpasses the mean \(FFT\) intensity plus \(K_{FFT}\) times its standard deviation \(\sigma\) (for example, for the 2-10 s band, the anomaly is defined if \(FFT(2s\leq p\leq 10s)>FFT+K_{FFT}\cdot\sigma(FFT)\). This criterion on \(K_{FFT}\) is in addition to Step 4 about \(k_{t}\), but it is not applied for the basic anomaly extraction (i.e., when no frequency analysis is made). We recognize that the residuals with respect to the spectral mean are not Gaussian, but it is feasible to introduce the standard deviation simply as an operational measure of the data dispersion around the mean of the spectrum. Even if other quantities (e.g., the absolute mean deviation or the largest deviation) could be used, we found that the results are practically the same, so we preferred to use the standard deviation that can be easily and quickly estimated. The two datasets (satellite anomalies and earthquakes) are investigated by the World-wide Statistical Correlation (WSC) algorithm fully described in [29] that is only briefly summarized here. The WSC algorithm first extracts the subset of anomalies and earthquakes to be correlated by a Superposed Epoch and Space Approach (SESA). The anomalies are selected by applying several thresholds. In particular, we only select tracks in geomagnetic quiet time (i.e., when the geomagnetic indices satisfy the following conditions: \(\mid\)Dst \(\mid\)\(\leq\) 20 nT and \(\text{a}_{\text{p}}\)\(\leq\) 10 nT), with all samples acquired in nominal satellite conditions (checked by the \"Quality Flags\" of the _Swarm_ satellite, according to the respective operating manuals), and put conditions over the rms of the window and, for some analyses, an additional threshold on the frequency. For each earthquake of magnitude \(M\), the anomalies inside its preparation area, as defined by a circular area with a Dobrovolsky's radius (in km: \(10^{(0.43\cdot M)}\)[36]) that occurred in the analyzed window (for example, from 500 days before the earthquake until 500 days after it), are extracted. If an anomaly is univocally associated with one earthquake, this couple (anomaly-earthquake) is defined.