## 2 Data and Methods
### Spectral Parameters Analysis  
The spectrum of a seismic signal can be parametrized by using the statistical moments of the power spectrum [20]. The statistical moments \(\lambda_{n}\) of the seismic power spectrum \(G\) were evaluated using the relation, Equation (1):  
\[\lambda_{n}=\int_{f_{1}}^{f_{2}}G\cdot f^{n}df \tag{1}\]  
where the index n corresponds to the n-th moment of the seismic power spectrum \(G\). The integration is evaluated in frequency domain in the range \(f_{1}\)-\(f_{2}\). For a single spectrum, the first three statistical moments \(\lambda_{0}\), \(\lambda_{1}\) e \(\lambda_{2}\) were used to evaluate the two spectral parameters,  
Figure 1: The seismic monitoring network operating on Mt. Vesuvius during the last years. The green triangles indicate the position of stand-alone seismic stations, while the blue symbols show the position of permanent network stations. The green symbols with a ticker border (BKWG, BKSG and SVAG) are the stations used in this work.  
i.e., central frequency \(\Omega\) and shape factor \(\delta\), following the relationships (Equations (2) and (3) respectively):  
\[\Omega=\sqrt{\frac{\lambda_{2}}{\lambda_{0}}} \tag{2}\]  
\[\delta=\sqrt{1-\frac{\lambda_{1}^{2}}{\lambda_{2}\lambda_{0}}} \tag{3}\]  
The central frequency \(\Omega\) is a measure of the frequency where the signal power is higher, while the shape factor \(\delta\) indicates the dispersion of the power spectral density around the central frequency [21]. The latter takes values between 0 and 1, with higher values corresponding to larger bandwidths. The power spectrum \(G\) and then the spectral parameters were calculated on adjacent 30 s signal windows for the continuous seismic signals. For any seismic station, the power spectrum parameters in Equations (2) and (3) were evaluated, respectively, for each component of ground motion, and then the average central frequency and the average shape factor among the three components were computed. In this way, for a given time window, we obtained a pair of values \(\Omega\) and \(\delta\) for each seismic station. A duration of 30 s corresponds approximately to an earthquake of M\({}_{\text{D}}\) = 1.3. We considered this window size a good compromise between the minimum detectable earthquake and the calculation cost of the procedure used. In particular, statistical moments were evaluated through Equation (1) in the frequency band 1-40 Hz. The lower integration limit was set at 1 Hz to reduce the effects of low-frequency seismic noise, which is often quite high at a frequency below 1-2 Hz in broadband signals, while the upper integration limit is given by the upper frequency limit of the flat bandwidth of the analyzed data. The calculation of the statistical moments was carried out on the continuous signals recorded at the BKSG, BKWG and SVAG sites for the years 2019, 2020 and part of 2021. In order to highlight potentially interesting signals for the search of LF events, signal windows were selected according to the results obtained by [PERSON] et al. [16]. Specifically, signal windows characterized by a central frequency \(\Omega\) and a shape factor \(\delta\) in the ranges 3 Hz \(\leq\Omega\leq\) 6 Hz and 0.45 \(\leq\delta\leq\) 0.65, respectively, were considered in the present study. With this procedure, we obtain independent values of \(\Omega\) and \(\delta\) for any analyzed seismic stations.