Application of space-time spectral analysis for detection of seismic waves in gravitational-wave interferometer
Robert Szymko, Mateusz Denys, Tomasz Bulik, Bartosz Idźkowski, Adam Kutynia, Krzysztof Nikliborc, Maciej Suchi?ski
AApplication of space-time spectral analysis for detection of seismic waves in gravitational-waveinterferometer
Robert Szymko , Mateusz Denys ∗2 , Tomasz Bulik , BartoszIdźkowski , Adam Kutynia , Krzysztof Nikliborc , and MaciejSuchiński Faculty of Physics, University of Warsaw, Pasteura 5, PL-02093 Warsaw, Poland AstroCeNT — Particle Astrophysics Science and Technology Centre InternationalResearch Agenda, Nicolaus Copernicus Astronomical Center, Polish Academy ofSciences, Rektorska 4, PL-00614 Warsaw, Poland Astronomical Observatory, University of Warsaw, Ujazdowskie Ave. 4, PL-00478Warsaw, Poland
February 9, 2021
Abstract
Mixed space-time spectral analysis was applied for the detection of seis-mic waves passing through the west-end building of the Virgo interferome-ter. The method enables detection of every single passing wave, includingits frequency, length, direction, and amplitude. A thorough analysis aimedto improving sensitivity of the Virgo detector was made for the data gath-ered by 38 seismic sensors, in the two-week measurement period, from 24January to 6 February 2018, and for frequency range 5–20 Hz. Two domi-nant seismic-wave frequencies were found: 5.5 Hz and 17.1 Hz. The possiblesources of these waves were identified, that is, the nearby industrial complexfor the frequency 5.5 Hz and a small object 100 m away from the west-endbuiding for 17.1 Hz. The obtained results are going to be used to providebetter estimation of the newtonian noise near the Virgo interferometer. ∗ e-mail: [email protected] a r X i v : . [ a s t r o - ph . I M ] F e b Introduction
First events detected by gravitational-wave (GW) detectors in the past decadegave us hope to answer the most fundamental questions of astronomy and physics,concerning, e.g., evolution of galactics and black holes, gravitational colapse ofsupernovae, or limits of general relativity [1–6]. However, noise reduction in thesedetectors is most complex task, due to very low amplitude of the incoming signals.Seismic characterization of the vicinity of the GW detectors is necessary for thenoise reduction, especially for filtering newtonian noise coming from fluctuatinggravitational forces caused by density perturbations in the Earth [7–9]
Space-time spectral analysis (or the frequency-wavenumber method) is a straight-forward extension of harmonic analysis to more than one dimension. Apart froma temporal dimension, additional spatial dimensions are introduced. It is mostsuitable if the spatial dimensions are cyclically continuous, or at least have fixedboundaries. In such cases we can look for modes of variability in which spatialscales have particular temporal scales. If the behavior is harmonic (wavelike),then we expect space-time spectral analysis to isolate any such modes that arepresent in the data. For instance, if we make space-time spectral analysis of astringed instrument, we expect to find a definite relationship between the lengthscales and the time scales of the oscillations [10].The method presented in this article is based on the paper [11] (and refs.therein), where seismic noise structure has been estimated using arrays of detec-tors (see also [12]). This method is one of the most popular methods of seismic-wave detection to the present day [13–15]. It is straightforward and provides alot of useful information about the seismic activity at the same time. Some morecomplex alternatives to this method are, e.g., High-Resolution method of Capon[16], the broad-band frequency-wavenumber spectral method of Nawab et al. [17],or the Multiple Signal Classification scheme [18] (see also [19, 20]).Space-time spectral analysis has been used before to obtaining the Rayleighwave phase velocity dispersion curves [21]. It is also a widely explored methodin the atmospheric sciences [22–26]. In [27, 28] space-time spectral analysishas been used for the description of zonal component of winds over the equator.Moreover, [29] applied space-time spectral analysis to modeling rainfall fields. Asfar as only spatial Fourier transform is concerned, Kim and Barros [30] used itstwo-dimensional version to analyse structure of soil moisture fields in differenttime moments, considered separately. That is, a version with one spatial dimension plus temporal dimension. Investigations of the seismic characteristics of the Virgo site have been madepreviously in [31, 32]. For some recent works about seismic-noise measurementsfor gravitational-wave detectors see, e.g., [33–36].This paper is structured as follows. The system of seismometers that we usedfor data aquisition is presented in Sec. 2. Then, we provide information aboutdata analysis process (Sec. 3). The results of this analysis are briefly described inSec. 4. Finally, some concluding remarks are given in Sec. 5.
During the study a system of seismic detectors in the west-end building (WEB)of the Virgo interferometer was examined. Location of the WEB as a part ofthe Virgo detector is shown on Fig. 1. The system consists of 38 seismometers, http://public.virgo-gw.eu/advanced-virgo/ (0 ,
4f the coordinate system shown on Fig. 2 is a beam spliter in the central building.The axes coincide with the arms of the detector, therefore, the axes form an angleca. ◦ with cardinal directions.The data that were used for the calculations were collected during a two-weekmeasurement period, from 24 January to 6 February 2018. The data for eachdetector are values of voltage changes due to its vertical motion, collected withthe sampling frequency 500 Hz, continuously. Thus, the system measures surfaceRayleigh waves. Not all (256) 1-hour data files generated in the measurementperiod are complete. Actually, the real measurement time is 250 h ( ∼ . days).During the experiment a series of each sensor’s 10 ’zero’ measurements one afteranother was assumed being empty-data period for this sensor. Collectively, thisis . of measurement time for all sensors. On Fig. 3 empty periods for allsensors are presented. The breaks in data collection for all sensors are not showed,therefore, hours on the horizontal axis do not correspond to real time.Although the sensors measure ground motions, actually the system was cre-ated in order to estimate newtonian noise in the Virgo measurement signal, i.e., anoise coming from fluctuating gravitational forces caused by density perturbationsin the Earth, which can be estimated using seismic-noise signal [7, 8]. Newtoniannoise has to be subtracted from the signal in any ground-based gravitational-wavedetector (cf. [34, 37]). The goal of this analysis is to improve sensitivity of theVirgo for low frequencies. For the above mentioned data, the following formula for the mixed space-timeFourier transform was applied,FT ( k x , k y , ω ) = (cid:90) t (cid:90) π − π (cid:90) π − π e − π i ( xk x + yk y + tω ) ψ ( x, y, t ) d x d y d t, (1)where ψ ( x, y, t ) is a measurement value of the seismometer placed in ( x, y ) intime moment t , ω is a frequency of the detected wave, and k x , k y are the respec-tive components of its two-dimensional wave vector k . Actually, a discretizedversion of Eq. (1) was used for ∆ t = 0 . s. The value of FT ( k x , k y , ω ) providesinformation about the amplitude of the particular wave. Notably, both spatialdimensions should not exceed ( − π, π ) range, so they should be normalized to thatrange (see below).In general, the space-time Fourier transform consists of wave vectors for eachwave detected in the signal ( k x , k y axes) together with the frequency of that Because the system of seismometers measure surface waves, z axis is not included in theabove formula. ω axis). Thanks to that we can obtain length, direction, frequency andamplitude of every wave that passed through the measuring system during theconsidered time period. Notably, Rayleigh surface waves are expected to prevailin the seismic noise in the vicinity of the Virgo building [38]. Therefore, we canneglect vertical differences between the spatial positions of the seismometers ( z axis).Remarkably, space-time Fourier transform, as any other multi-dimensionalFourier transform, can be calculated in an independent way, cf. Eq. (1). Thatis to say, having the input data of ( x, y, t ) -points, we can calculate, firstly, onlythe temporal part [obtaining ( x, y, ω ) -points] and secondly, using this result, wecan subsequently calculate the spatial part, obtaining the final result, ( k x , k y , ω ) -points. We may also make that inversely, i.e., calculating first the spatial partand then the temporal part to obtain the final result.In the approach presented herein, initially, the average value of the signalwas subtracted from the input, for each detector separately. After that, thewhole signal was divided into 20-second pieces (time frames), each of them with N = 10 000 samples. Thanks to dividing the data into short pieces, we candetect short wave impulses, lasting only a couple of seconds. Due to the fact thatthe space-time method, as any other Fourier-based method, uses finite sets ofdata, the questions of spectral leakage and, consequently, the necessity of applyingappropriate window functions should be considered, for some artificial componentsmay appear in the final result. For that purpose the Hann window function wasused, w ( n ) = 12 (cid:20) − cos (cid:18) πnN − (cid:19)(cid:21) , (2)where n indexes each measurement in the N -sample time frame. The seismome-ters are 24-bit and measure voltage changes in the range between − V and 5 Vwith sensitivity 77.3 V/m/s. Therefore, the input signal of the seismometerswas converted from volts to m/s, by using the appropriate conversion coefficient / (2 · . .After the initial preparation of data, explained above, the time-frequency part,for every time frame and for each detector separately, was calculated, using FFT(Fast Fourier Transform) algorithm. After that, the result was multiplied by theappropriate detector response function, which depends on the frequency value ω .Subsequently, the spatial part was calculated, using the values obtained in thefirst step. To calculate the spatial part in the second step, variables ( x, y ) were By calculating an inverse Fourier transform of each ASD peak in the final result. However,we have to consider in this case all rescalings done during the process of the FT calculation —see below in this section. u = s f x and v = s f y , using the same scaling factor, s f = π max {| x | max , | y | max } , (3)so that both spatial dimensions were normalized to ( − π, π ) range without dis-torting the waves. The ( − π, π ) range of the spatial data was required by theNUFFT (Non-Uniform FFT) algorithm that was used subsequently. A version ofthe NUFFT algorithm for two dimensions was applied [39]. NUFFT is the versionof FFT algorithm for nonuniform samples, as the locations of the seismometerswere not uniform. Finally, after calculating the spatial part, the variables werechanged back to k x = s f k u , k y = s f k v , and three-dimensional complex modulus ofthe spectrum was calculated to obtain the 3D amplitude spectral density (ASD):ASD ( k x , k y , ω ) = | FT ( k x , k y , ω ) | , (4)where FT ( k x , k y , ω ) is given by Eq. (1). In order to test the algorithm presented above, the following simulation was run.A simulated set of data (adopting the same ( x, y ) positions of the seismometers)with the superposition of two sinusoidal waves: ψ ( x, y, t ) = A cos( k x x − ω t ) + B cos( k y y + ω t ) (5)is created, where k x = 0 . m − , k y = 0 . m − , ω = 3(2 π ) Hz, ω = 5(2 π ) Hz,and amplitudes of both waves A = B = 1000 V. Figure 4 shows the obtained space-time amplitude spectral density for theseexemplary data, in a form of colormap, e.g., for 3, 4, and 5 Hz. Additionally,a rank of ASD values is shown for frequency 3 Hz. As expected, the space-timeamplitude spectral density for such strong signal consists of two outstanding peaksfor frequency 3 Hz and 5 Hz [Figs. 4(a) and 4(c)], corresponding to both partsof the superpositioned wave, with proper wave-vector and frequency values. Thepeak height for 5 Hz is about two times smaller than for 3 Hz, due to applyingdetector response function, which partially dampens higher frequencies (cf. Sec.3). However, both peaks are outstanding compared to the other spatial modes forfrequencies 3 and 5 Hz, respectively: they are more than 10 times larger than thestandard deviation for the whole plot for the particular frequency [on the rankplot for 3 Hz, Fig. 4(a) on the right, the outstanding k -value is denoted by the red To mimic the seismometer measurements in volts, approximated to integer values, due tothe fact that the real seismometer measurements in volts were integer values. k =0 . m − , k = 0 . m − , ω = 3(2 π ) Hz, ω = 5(2 π ) Hz, A = B = 1000 V) for(a) 3 Hz (on the left; on the right, additionally, a rank of ASD values in semi-logscale is shown, where the outstanding point corresponding to the yellow peak onthe left is marked in red), (b) 4 Hz, (c) 5 Hz.9ot in the bottom right corner]. This is also true for the neighboring frequencyvalues, since both peaks have a wide profile in the frequency domain — peaksare visible for a large frequency range around 3 and 5 Hz, respectively, with thevalue of the peak larger than 10 standard deviations of all spatial modes in theplot. However, the result for, e.g., 4 Hz is much weaker than for 3 and 5 Hz, andcompared with them shows no signal at all [Fig. 4(b)].Therefore the importance of the ranking-plot representation shown on Fig. 4(a)is that it allows the distinction of the essential points in the amplitude spectraldensity plot. Consequently, the future signal processing in order to separate thereal seismic-wave signal from the noise in the spectral-density data may use theranking-plot analyses. Nevertheless, for the real data we may expect the obtainedpeaks of ASD being weaker than in the above test example, due to the fact thatthe test data contained strong input signal without any additional noise, whereasfor the real data the input signal is expected to contain some noise.
After the successful test of the algorithm on the artificial data, the analysis ofthe real empirical data collected from the Virgo West-End building was made.In the initial step, the seismometer indications in volts were converted to theappropriate values in m/s. Then, the three-dimensional transforms for these datawere calculated using FFT and NUFFT algorithms (cf. Sec. 3). Finally, the space-time amplitude spectral density was calculated for each 20-second period and eachfrequency. An example plot is usually characterized by one main maximum — cf.Fig. 5.The next step was finding a position ( k x , k y ) of maximal ASD value in time,for each frequency from the range 5 to 20 Hz, step 0.1 Hz. Examples of pre-liminary results for 20-second pieces of seismic data, taken from the same hourof measurement, are shown in Fig. 6. There are dominant waves for particularfrequencies, as expected. On first plot for 5 Hz [Fig. 6(a)] there is a strong max-imum (north-east wave direction). On Fig. 6(b) for 10 Hz there is much weakerand fuzzy maximum. For (c) and (d) (12 Hz and 15 Hz, respectively) there arenot any signals visible, but rather some artifacts present due to inaccuracy of themeasurement or the NUFFT method used.On Fig. 7 ASD values for ( k x , k y ) = (0 , and for each frequency, averagedfor one example hour of measurement, are shown. Observed maxima of ASD arevery sharp, therefore the step 0.1 Hz is taken for this plot.10igure 5: ASD value in ( k x , k y ) space given by Eq. (1) for frequency 5.5 Hz, 10thminute of the first hour of measurements (marked also using colorscale on theleft). 11igure 6: Colormap of ASD in ( k x , k y ) space for 20-second pieces of seismic datafrom 2018-01-24, 13:30, and for: (a) 5 Hz, (b) 10 Hz, (c) 12 Hz, (d) 15 Hz.12igure 7: ASD for ( k x , k y ) = (0 , and frequencies from 5 to 20 Hz, averaged overan example hour of measurement. Resolution of wave-vector space taken for our calculation is δ k = 0 . m − , thus,for each frequency 180 matrices × were obtained. In order to determinedominant waves the following function (paraboloid) was fitted to the space-timeFourier spectra: fit ( u, v ) = A − uσ u − vσ v , (6)where A is maximal amplitude for a given matrix × , σ u and σ v are standarddeviations of u and v , respectively, u = ( k x − k x ) · cos( α ) − ( k y − k y ) · sin( α ) , (7) v = ( k x − k x ) · sin( α ) + ( k y − k y ) · cos( α ) , (8) k x , k y are coordinates of maximal ASD value, α is inclination of the paraboloidaxis with respect to the X axis, and k x , k y are the matrix coordinates.In order to increase the resolution in wave-vector space, for each matrix anarea of × points close to the maximum was inspected, with an accuracy of max { ASD ( k x , k y ) } increased by two orders of magnitude. Then, the logarithm, L , of the likelihood function: L = (cid:88) i,j =1 [ fit ( k x , k y , k x , k y , σ u , σ v , α, A ) − λ ( i, j ) /σ L ] (9)13as calculated, where λ ( i, j ) = ASD ( k x , k y , ω ) is a value of ASD-matrix elementfor a particular ω , and σ L = 0 . · A is assumed uncertainty of the experiment. The values of σ u and σ v in Eq. 6 were checked in the range 50 to 150 m. Afternormalization, collected values of L were aggregated in the following way L x,y ( k x , k y ) = (cid:88) σ u ,σ v ,α L ( k x , k y , σ u , σ v , α ) (10) L x ( k x ) = (cid:88) k x L ( k x , k y ) (11) L y ( k y ) = (cid:88) k y L ( k x , k y ) (12)Finally, the cumulative distribution function and the coordinates of ASD max-imum for each frequency and each 20-second period of the measurement wereobtained. In this section the final results for all data gathered between 24-01-2018 and 06-02-2018 are presented. Fig. 8 shows dominant frequencies for the whole period.Apparently, there are two main frequencies for ASD: ca. 5-6 Hz and ca. 17 Hz.The first one is characterized by strong day-night periodicity. On Figs. 9 and10 histograms of the data are shown. Fig. 9 shows a distribution of max(ASD)values. Fig. 10 shows a distribution of dominant frequencies. Two main dominantfrequencies are 5.5 Hz and 17.1 Hz.In the next step, a detailed analysis was made in order to characterize seismicwaves in the vicinity of the Virgo WEB. Using methodology from Sec. 3.3 exactcoordinates of max(ASD) in the wave-vector space, for particular frequencies and20-second periods of the measurement, were found. Collected results are shownon Fig. 11 for 17.1 Hz and Fig. 12 for 5.5 Hz. For the sake of clarity, the resultsare projected on XZ , Y Z , and XY surfaces, where X and Y reflect k x and k y ,respectively, and Z reflects ASD. Due to presence of noises in the measurementsystem and inaccuracy of our method, a lower limit for max(ASD) was set to . of the maximal ASD value for all data, which was equal to . m / s / .Figs. 13 and 14 for frequencies 17.1 Hz and 5.5 Hz, respectively, are histogramsof wave direction (upper parts) and wavenumber (lower parts). o on upper plots A change of σ L does not change coordinates of the ASD maximum, just its value. RHS of Eq. (10) does not include k x , k y as they reflect indices of summation i, j in Eq. (9),and A , as it was constant for all calculated L values. Counts for zero value on both plots are a consequence of quality of the data and should notbe considered at all — cf. Fig. 3. · − ).15igure 10: Histogram of dominant frequencies (bin width 0.13).does not indicate cardinal direction, but direction of X axis in the seismometercoordination system shown on Fig. 2. For frequency 17.1 Hz results indicateclearly the direction − ± . o and wavenumber | K . | = 0 . . These valuescorrespond to the seismic wave arriving from northeast with velocity ± m/s. Results for 5.5 Hz are more ambiguous and finding a strict wave directionin this case is impossible. The most frequent direction is ca. − ± . o withwavenumber | K . | = 0 . . These values correspond to the seismic wave arrivingfrom northwest with velocity ± m/s.On Fig. 15 the map of the vicinity of the Virgo detector with directions ofwaves propagation is shown. Red lines is an assumed coordination system, shifted . o clockwise to true north. Blue lines indicate direction of 5.5-Hz wave witherrors. The area in blue rectangle, i.e., an industrial complex 2.5 km away fromthe WEB, may be a source of this wave. Black color indicate direction of 17.1-Hzwave. In top right corner of Fig. 15 a zoom of the vicinity of the WEB is shown,with a potential source of 17.1-Hz wave marked (black rectangle), ca. 100 m awayfor the WEB. Using the data from the experiment conducted between 24-01-2018 and 06-02-2018 collected by the system of seismometers inside the west-end building of the16igure 11: Locations of max(ASD) in wave-vector space for frequency 17.1 Hzprojected onto (a) XZ , (b) Y Z , (c) XY surface where X and Y reflect k x and k y , respectively, and Z reflects ASD. Colors indicate ASD values.17igure 12: Locations of max(ASD) in wave-vector space for frequency 5.5 Hzprojected onto (a) XZ , (b) Y Z , (c) XY surface where X and Y reflect k x and k y , respectively, and Z reflects ASD. Colors indicate ASD values.18igure 13: Histograms of direction (top) and module (bottom) of the wave vectorfor frequency 17.1 Hz. Bin width: . o (top), · − . o (top), · − × km). Red lines indicatecoordination system consistent with that shown on Fig. 2. The directions of wavespropagation for 5.5-Hz wave with error ± . o (blue arrow + side lines) and for17.1-Hz wave with error ± . o (black arrow + side lines). Rectangles indicate thepossible wave sources. The area close to the 17.1-Hz wave source is zoomed andshown in the top right corner of the figure.20irgo gravitational interferometer, we made a seismic analysis of the nearby area.The seismic noise (currently compensated) corresponds to the changes in thegravitational field that can be detected by the Virgo. These fluctuations shouldbe also identified in order to increase sensitivity of the GW measurement system.Particularly, for a frequency below 20-30 Hz they are recognized as newtoniannoise [37, 40, 41].In this work, we searched for the seismic waves with a frequency from 5 to20 Hz. Two dominant frequencies were detected and identified with two mainsources, the first one emitting a continuous wave, and the second one emittingwave cyclically. First observed dominant frequency is 17.1 Hz. The directionof this wave is − ± . o , i.e., it was arriving from northeast, with velocity ± m/s. On Fig. 15 in the top right corner a small object, ca. 100 m awayfrom the WEB is shown. It is located exactly at the found wave direction andcan be a source of 17.1-Hz wave.However, the source of the 17.1-Hz wave can be also internal, located inside thebuilding. Therefore, another possible interpretation of this wave is a resonanceof the whole building, which is excited by vacuum pumps, cooling fans, or otherdevices [36].Seismic waves at the frequency 5.5 Hz coincide with daily human activity.They dominate between 6 a.m. and 6 p.m. every day and move with velocity ± m/s. A direction of this waves is not strict, ca. − ± . o . This iscaused by a supposed origin of the waves, which is human activity in the nearby in-dustrial complex. The study made before, indicated another antropogenic sourceof the noise at the frequency 4-10 Hz: three nearby bridges [42]. The velocitiesgiven above are typical for dispersion media, i.e., the wave with a lower frequencymoves faster, and agree with the predictions of one of the models from [8]. Thetypical velocity of Rayleigh waves in metals is 2–5 km/s, while in soil 50–700 m/sfor shallow waves (under 100 m) and 1.5–4 km/s for depth over 1 km [43]. The method of seismic-wave detection in the Virgo WEB presented above, is ableto detect particular waves that pass through the building. The main advantages ofthe method is that it gives the full information about the waves that come throughthe measuring system (frequency, length, direction, and amplitude) with low lossof the information from the seismometers comparing, e.g., to the method usingcoherence values between each pair of the detectors during a certain measurementperiod, for example, one hour [36, 44, 45]. Our method is more sensitive andefficient, and it is most useful for the characterization of seismic field, as it wasshown in this work. 21he idea of minimization of newtonian noise in the Virgo GW signal is basedon monitoring sources of gravitational field perturbations. These perturbationsare caused by local changes of density in the vicinity of the mirrors. A commonmethod of filtering newtonian noise is based on Wiener filters, which have beenalready applied to reduction of other noises in GW detectors. The parametersof seismic waves found in this work can be used to fit correlation models of theWiener filter [46].As the system of mirrors in the Virgo is complex (see Fig. 1), the analysisconducted herein shall also be made for other parts of the detector. For the west-end building itself the data were collected for more than a year; in the north-endbuilding a system of sensors is going to be installed as well. Complementarily tothe system of seismometers, infrasound microphones are going to be used for aprecise determination of the sources of gravitational field fluctuations caused bypressure changes. Particularly, a detailed analysis of the frequencies below 10 Hzwould be vital, as for 1.7 Hz a significant noise coming from the nearby wind farmis observed [32].One of the problems that needs solving when we use the space-time spectralanalysis method is how to distinguish between the real waves and the noise. Forinstance, some threshold values may be used for this purpose [47, 48]. However,in this article we focused on the dominant frequencies observed in the obtainedsignal. Moreover, some problems with the spatial aliasing in the results should beresolved, e.g., by applying the window function also for the spatial dimensions (cf.Sec. 3). On top of that, we may consider the use of overlapping time periods whencalculating the spectra, to improve the accuracy and relevance of the results. Thefurther studies will contribute to the improvement of quality of the informationabout the gravitational-wave sources like binary black holes and neutron stars.
Acknowledgments
M.D. is grateful for inspiring discussion with Krzysztof Lorek and Jan Harms.The authors are grateful to colleagues from the Virgo collaboration and the VirgoNewtonian noise team for providing the seismic array data and for fruitful dis-cussions. The work was supported by the TEAM grant from the Foundation forPolish Science TEAM/2016-3/19 and by the grant
AstroCeNT — Particle As-trophysics Science and Technology Centre International Research Agenda carriedout within the International Research Agendas programme of the Foundation forPolish Science financed by the European Union under the European RegionalDevelopment Fund. 22 ppendix A Coordinates of the seismic sensors
No. x [ m ] y [ m ] z [ m ] No. x [ m ] y [ m ] z [ m ] Table 1: Coordinates of the seismometers with respect to the Virgo test mass.
References [1] B.P. Abbott et al. , Phys. Rev. Lett. , 6 (2016).[2] A. Sesana,
Phys. Rev. Lett. , 23 (2016).[3] S. Endlich et al. , J. High Energy Phys. , 9 (2017).[4] M. Sasaki et al. , Class. Quantum Gravity , 6 (2018).[5] M. Miller, M. Coleman, and N. Yunes, Nature , 7753 (2019).[6] B.P. Abbott et al. , Phys. Rev. D , 10 (2019).[7] P. R. Saulson,
Phys. Rev. D , 4 (1984).238] S. A. Hughes and K. S. Thorne, Phys. Rev. D , 12 (1998).[9] M. Coughlin et al. , Class. Quantum Gravity , 24 (2016).[10] D. L. Hartmann, ATM S, (2008).[11] R. T. Lacoss, E. J. Kelly, and M. N. Toksöz, Geophysics , 1 (1969).[12] J. Gazdag and P. Sguazzero, Proceedings of the IEEE , 10 (1984).[13] H. Sato, M. C. Fehler, and T. Maeda, Seismic wave propagation and scatter-ing in the heterogeneous earth (Springer Science & Business Media, 2012).[14] S. Foti et al. , Surface wave methods for near-surface site characterization (CRC press, 2014).[15] S. Foti et al. , Surv. Geophys. , 6 (2011).[16] J. Capon, Proceedings of the IEEE , 8 (1969).[17] S. Nawab, F. Dowla, and R. Lacoss, IEEE Trans. Acous. Speech Sig. Proc. , 5 (1985).[18] R. Schmidt, IEEE Trans. Antennas Propagation , 3 (1986).[19] S. Rost and C. Thomas, Rev. Geophys. , 3 (2002).[20] Observatory Seismology: A Centennial Symposium for the Berkeley Seis-mographic Stations , edited by J. J. Litehiser (University of California Press,2018).[21] G. Mainsant et al. , J. Geophys. Res. Earth Surf. , F1 (2012).[22] S. K. Kao,
J. Atmos. Sci. , 1 (1968).[23] Y. Hayashi, J. Meteorol. Soc. Jpn. Ser. II
1, (1982).[24] M. L. Salby,
J. Atmos. Sci. , 11 (1982).[25] D. L. Wu, P. B. Hays, and W. R. Skinner, J. Atmos. Sci. , 20 (1995).[26] H. H. Hendon and M. C. Wheeler, J. Atmos. Sci. , 9 (2008).[27] Y. Hayashi, J. Meteorol. Soc. Jpn. Ser. II , 4 (1977).[28] Y. Hayashi, J. Atmos. Sci. , 6 (1979).2429] C. De Michele and P. Bernardara, Atmospheric Res. , 1-4 (2005).[30] G. Kim and A. P. Barros, Remote Sens. Environ. , 2-3 (2002).[31] F. Acernese et al. , Class. Quantum Gravity , 5 (2004).[32] G. Saccorotti, D. Piccinini, L. Cauchie, and I. Fiori, Bull. Seismol. Soc. Am. , 2 (2011).[33] L. Somlai et al. , Acta Geod. Geophys. , 2 (2019).[34] P. Ván et al. , Eur. Phys. J. ST , 8 (2019).[35] N. Mukund et al. , Class. Quantum Gravity , 8 (2019).[36] M. C. Tringali et al. , Class. Quantum Gravity , 2 (2019).[37] J. Harms and H. J. Paik, Phys. Rev. D , 2 (2015).[38] A. Singha, S. Hild, and J. Harms, Class. Quantum Gravity , 10 (2020).[39] A. Dutt and V. Rokhlin, SIAM J. Sci. Comput. , 6 (1993).[40] J. Harms, Living Rev. Relativ. , 6 (2019).[41] J. C. Driggers, J. Harms, and R. X. Adhikari, Rana X., Phys. Rev. D , 10(2012).[42] F. Acernese et al. , Class. Quantum Gravity , 5 (2004).[43] W. M. Telford, L. P. Geldart, and R. E. Sheriff, Applied geophysics (Cam-bridge University Press, 1990), p. 149.[44] R. S. Harichandran and E. H. Vanmarcke,
J. Eng. Mech. , 2 (1986).[45] M. Coughlin et al. , J. Geophys. Res.-Sol. Ea. , 3 (2019).[46] M. Coughlin et al. , Class. Quantum Gravity , 24 (2016).[47] M. Denys, T. Gubiec, R. Kutner, M. Jagielski, and H. E. Stanley, Phys. Rev.E , 4 (2016).[48] J. Schäfer, J. Atmos. Sci.36