# Extension of the spatial autocorrelation (SPAC) method to mixed-component correlations of surface waves,

M. Hanay, T.D. Mikesell, K. van Wijk, H. Nakahara

► To cite this version:

M. Hanay, T.D. Mikesell, K. van Wijk, H. Nakahara. Extension of the spatial autocorrelation (SPAC) method to mixed-component correlations of surface waves. Geophysical Journal International, 2012, 191 (1), pp.189-206. 10.1111/j.1365-246X.2012.05597.x. http://hal-00812671

HAL Id: hal-00812671

https://hal.science/hal-00812671

Submitted on 17 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Extension of the spatial autocorrelation (SPAC) method to mixed-component correlations of surface waves

Matthew M. Haney,<fnref n="1" /> T. Dylan Mikesell,<fnref n="2" /> Kasper van Wijk<fnref n="3" /> and Hisashi Nakahara<fnref n="4" />

Accepted 2012 June 28. Received 2012 June 19; in original form 2011 November 3

## SUMMARY

Using ambient seismic noise for imaging subsurface structure dates back to the development of the spatial autocorrelation (SPAC) method in the 1950s. We present a theoretical analysis of the SPAC method for multicomponent recordings of surface waves to determine the complete 3 × 3 matrix of correlations between all pairs of three-component motions, called the correlation matrix. In the case of isotropic incidence, when either Rayleigh or Love waves arrive from all directions with equal power, the only non-zero off-diagonal terms in the matrix are the vertical–radial (ZR) and radial–vertical (RZ) correlations in the presence of Rayleigh waves. Such combinations were not considered in the development of the SPAC method. The method originally addressed the vertical–vertical (ZZ), RR and TT correlations, hence the name spatial autocorrelation. The theoretical expressions we derive for the ZR and RZ correlations offer additional ways to measure Rayleigh wave dispersion within the SPAC framework.

Expanding on the results for isotropic incidence, we derive the complete correlation matrix in the case of generally anisotropic incidence. We show that the ZR and RZ correlations have advantageous properties in the presence of an out-of-plane directional wavefield compared to ZZ and RR correlations. We apply the results for mixed-component correlations to a data set from Akutan Volcano, Alaska and find consistent estimates of Rayleigh wave phase velocity from ZR compared to ZZ correlations. This work together with the recently discovered connections between the SPAC method and time-domain correlations of ambient noise provide further insights into the retrieval of surface wave Green's functions from seismic noise.

**Key words:** Time series analysis; Interferometry; Surface waves and free oscillations; Theoretical seismology; Wave propagation.

## 1 INTRODUCTION

Seismologists often encounter signals that do not exhibit distinct *P*- or *S*-wave phase arrivals. Examples of these signals include oceanic microseisms (Tanimoto & Alvizuri 2006; Muyzert 2007a,b), the Earth's hum (Rhie & Romanowicz 2004), glacial earthquakes (Ekström et al. 2003), volcanic tremor (Ferrazzini et al. 1991; Haney 2010), volcanic explosions (Chouet et al. 2005) and debris flows such as lahars (Zobin et al. 2009). For these cases, the seismologist's standard toolbox, based on picking the arrival times and amplitudes of body-wave phase arrivals, is not applicable. At best, for short duration signals such as regional earthquakes and teleseisms, passing surface waves can be identified and used to infer structure of the crust and upper mantle (Dorman & Ewing 1962). Recent advances in seismology, including ambient noise tomography (ANT), have pioneered the use of extended signals lacking body waves for the purpose of imaging Earth structure. ANT achieves this goal by correlating long-time recordings of ambient seismic noise within the frequency band dominated by the oceanic microseism (<1 Hz). Although the proliferation of ambient noise techniques has mostly occurred during the past decade, the idea of correlating long duration seismic signals from spatially extended sources goes back much further. In fact, the basic approach to interstation correlation originates over 50 yr ago, in seminal papers by Aki (1957, 1965). Aki's spatial autocorrelation (SPAC) method has subsequently been used to study shallow subsurface structure using wavefields composed of microseisms (Okada 2003) and volcanic tremor (Ferrazzini et al. 1991; Chouet 1996; Chouet et al. 1998; Saccorotti et al. 2003). Okada (2003) provided a comprehensive review of the SPAC method, also known as the microtremor method.

<footnote n="1">¹Alaska Volcano Observatory, U.S. Geological Survey, Anchorage, AK 99508, USA. E-mail: mhaney@usgs.gov</footnote>

<footnote n="2">²Géoazur, Centre National de la Recherche Scientifique (UMR 6526) Observatoire de la Côte d'Azur, Université de Nice Sophia-Antipolis, Valbonne, France</footnote>

<footnote n="3">³Department of Geosciences, Boise State University, Boise, ID 83725, USA</footnote>

<footnote n="4">⁴Department of Geophysics, Graduate School of Science, Tohoku University, Aoba-ku, Sendai 980-8578, Japan</footnote>
Recently, Aki's work has been revisited in light of advances made in the use of ambient noise following the groundbreaking work of Campillo & Paul (2003). Connections between time-domain ambient noise correlations and the SPAC method have been explored by several authors (Chávez-García & Luzón 2005; Chávez-García et al. 2005; Nakahara 2006; Ekström et al. 2009; Prieto et al. 2009; Tsai 2009; Harmon et al. 2010; Tsai & Moschetti 2010). In the presence of an isotropic noise field, Nakahara (2006) connected the Fourier transform of the time-domain ambient noise correlation and the spectral correlation coefficient from the SPAC method. Building on this connection, Ekström et al. (2009) showed that the zeros in the real part of the spectrum of averaged time-domain correlations could be associated with the spectral zeros predicted by Aki's SPAC theory. In the case of an anisotropic noise field, Nakahara (2006) derived expressions for time-domain ambient noise correlations drawing on results derived in the field of acoustics (Cox 1973; Teal et al. 2002). Through these studies, the SPAC method and time-domain correlations have been shown to be different but largely equivalent ways of using ambient noise to gain information on the structure of the subsurface.

Despite the connections forged between the two methods, comparisons between SPAC and time-domain correlations have thus far been limited to discussions of correlations between vertical–vertical (ZZ) components, with some discussion of radial–radial (RR) and transverse–transverse (TT) components (Chávez-García & Luzón 2005). Cho & Tada (2006) and Tada & Cho (2007) have extended the SPAC method under general conditions for three-component recordings. However, mixed components of polarized waves, such as ZR and RZ, have not been considered before. This is due in part to the fact that Aki (1957) only discussed ZZ, RR and TT correlations in the original paper on the SPAC method. In fact, the SPAC method has its name in part because it considers correlations between identical components, and therefore, a comparison between the complete Green's tensor of time-domain correlations and the SPAC method cannot be made. Even so, several authors have presented the full matrix of possible time-domain ambient noise correlations from three-component instruments (Campillo & Paul 2003; Stehly et al. 2007; Roux 2009; Durand et al. 2011; Roux et al. 2011). Recent work by Van Wijk et al. (2011) on the mixed components in the time domain suggests these independent results provide additional information that may be more robust than the diagonal terms. In this work, we extend the SPAC method to the complete matrix of correlation coefficients and compare the result to time-domain correlations between three-component instruments. Furthermore, we quantify the sensitivity of the different components of the correlation matrix to out-of-plane, directional ambient noise. To do so, we proceed in the way shown by Nakahara (2006) for unpolarized waves recorded on a single vertical component.

For isotropic incidence, our results for the off-diagonal components of the Rayleigh $\phi^R$ and Love $\phi^L$ wave correlation coefficient matrices, combined with Aki's results for the main diagonal, can be summarized as follows:

$$ \phi^R(r, \omega) = \begin{bmatrix} \phi_{ZZ}^R & \phi_{ZR}^R & \phi_{ZT}^R \\ \phi_{RZ}^R & \phi_{RR}^R & \phi_{RT}^R \\ \phi_{IT}^R & \phi_{TR}^R & \phi_{TT}^R \end{bmatrix} = P^R(\omega) \times \begin{bmatrix} J_0(\omega r/c_R) & -R J_1(|\omega| r/c_R) & 0 \\ R J_1(|\omega| r/c_R) & R^2 [J_0(\omega r/c_R) - J_2(\omega r/c_R)]/2 & 0 \\ 0 & 0 & R^2 [J_0(\omega r/c_R) + J_2(\omega r/c_R)]/2 \end{bmatrix} \quad (1) $$

$$ \phi^L(r, \omega) = \begin{bmatrix} \phi_{ZZ}^L & \phi_{ZR}^L & \phi_{ZT}^L \\ \phi_{RZ}^L & \phi_{RR}^L & \phi_{RT}^L \\ \phi_{IT}^L & \phi_{TR}^L & \phi_{TT}^L \end{bmatrix} = P^L(\omega) \times \begin{bmatrix} 0 & 0 & 0 \\ 0 & [J_0(\omega r/c_L) + J_2(\omega r/c_L)]/2 & 0 \\ 0 & 0 & [J_0(\omega r/c_L) - J_2(\omega r/c_L)]/2 \end{bmatrix}, \quad (2) $$

where $\omega$ is angular frequency, $r$ is radial distance, $c_R$ is Rayleigh wave velocity, $c_L$ is Love wave velocity, $R$ is the ratio of the horizontal-to-vertical motion of the Rayleigh waves, $P^R$ is the power spectrum of the Rayleigh waves, $P^L$ is the power spectrum of the Love waves and $J_0$, $J_1$ and $J_2$ are Bessel functions of the zeroth, first and second orders, respectively. Note that the correlation matrix for the Love waves lacks non-zero off-diagonal components in this case. Thus, Aki's theory for Love waves is complete for isotropic incidence. Moving beyond an isotropic wavefield, we also derive the Rayleigh and Love wave correlation coefficient matrices for generally anisotropic incidence (Appendix A). For anisotropic incidence, all off-diagonal components are non-zero in general for both Rayleigh and Love waves.

In the following, we present the theory for Rayleigh wave propagation in 1-D and 2-D and make connections with the time-domain Green's functions. We show the robustness of ZR correlations in the presence of out-of-plane ambient noise. With data from two broadband seismometers at Akutan Volcano, Alaska, we use ZR correlations to obtain Rayleigh wave dispersion curve measurements within the SPAC methodology that are complementary to measurements from ZZ correlations.

## 2 THE 1-D CASE FOR RAYLEIGH WAVES

We begin by analysing Rayleigh waves in a 1-D setting. Although the results we derive are not the ultimate findings of our study, the 1-D case has many conceptual similarities to the 2-D case. In this way, we can demonstrate the extensions we propose to the original SPAC method of Aki (1957) without some of the complexities of the 2-D case. Our derivation follows closely the development of SPAC in Aki (1957) and we adopt much of the notation used in that paper. We refer to the original manuscript for detailed discussions on the assumptions inherent to the method.
The 1-D case generally consists of Rayleigh waves propagating to the left or right within a single vertical plane. The horizontal $u_x$ and vertical $u_z$ displacements associated with the Rayleigh wave satisfy

$$u_z(x, t) = \sum_{n=-\infty}^{\infty} A_n \exp(i\rho_n x) \cos(c\rho_n t) + \sum_{n=-\infty}^{\infty} \frac{B_n}{c\rho_n} \exp(i\rho_n x) \sin(c\rho_n t) \quad (3)$$

$$u_x(x, t) = \sum_{n=-\infty}^{\infty} -RA_n \operatorname{sgn}(\rho_n) \exp(i\rho_n x) \sin(c\rho_n t) + \sum_{n=-\infty}^{\infty} R \frac{B_n}{c\rho_n} \operatorname{sgn}(\rho_n) \exp(i\rho_n x) \cos(c\rho_n t) \quad (4)$$

where c is the wave speed, R is the ratio of the horizontal-to-vertical motion and $\rho_n = 2\pi n/L$ is the wavenumber associated with the nth lateral mode over a finite interval of length L. Eq. (3) has the same form as the solution to the 1-D wave equation presented by Aki (1957). To analyse Rayleigh waves, we have added the coupled equation for $u_x$ in eq. (4). The Rayleigh waves considered here are assumed to consist of a single vertical mode, for example, the fundamental mode as in Aki (1957). Note that $u_x$ is different from $u_z$ in the following ways:

(i) It generally has a different amplitude, represented by the ellipticity R.

(ii) It is phase-shifted by 90° in time.

(iii) The sense of its motion changes whether the Rayleigh wave is propagating to the left or right.

Since $u_x$ and $u_z$ are real, the factors $A_n$ and $B_n$ are the complex conjugates of $A_{-n}$ and $B_{-n}$. Furthermore, as shown in eq. (4) of Aki (1957), the amplitudes of the different modes are assumed to be uncorrelated random variables (Lobkis & Weaver 2001):

$$\overline{A_n A_m} = |A_n|^2 \delta_{n,-m} \quad (5)$$

with an analogous relation for $B_n$. The symbol $\delta$ is the Kronecker delta function and the horizontal bar denotes averaging. Since they are uncorrelated random variables, the different lateral modes are taken to be in a state of white noise.

We first consider the spatial correlation of the vertical components, which was previously analysed by Aki (1957):

$$\phi_{ZZ}(\xi, t) = \overline{u_z(x, t) u_z(x + \xi, t)}. \quad (6)$$

Based on eq. (5) above, the spatial correlation is given by

$$\phi_{ZZ}(\xi, t) = \sum_{n=-\infty}^{\infty} \left[ |A_n|^2 \cos^2(c\rho_n t) + \frac{|B_n|^2}{c^2\rho_n^2} \sin^2(c\rho_n t) \right] \exp(i\rho_n \xi), \quad (7)$$

where we have used the property

$$\overline{A_n B_m} = -\overline{A_m B_n} \delta_{n,-m} \quad (8)$$

to cancel the two terms with dependence on time as $\sin(c\rho_n t)\cos(c\rho_n t)$. Note that this is different in form than the analogous property invoked in eq. (6) of Aki (1957), which was simply $\overline{A_n B_m} = 0$. In subsequent uses of these quantities, we omit the horizontal bar but the averaging is assumed. Thus, $\overline{Re(A_n B_{-n})} = \overline{Re(A_{-n} B_n)} = 0$ as in Aki (1957) but the imaginary parts of these terms are non-zero and opposite in sign.

Aki (1957) explains that the condition $\overline{A_n B_m} = 0$ physically means the initial distributions of displacement and particle velocity are independent. It is worth noting that the alternative condition we adopt in this paper, eq. (8), also satisfies this physical requirement. The condition in eq. (8) is less restrictive than the condition from Aki (1957): it implies that the sum of the terms $A_n B_{-n}$ and $A_{-n} B_n$ is zero instead of each term being equal to zero individually. This distinction was not important for Aki (1957) since mixed-components were not analyzed, but will play a major role here.

As pointed out by Aki (1957), the time dependence in eq. (7) can be removed if we assume that

$$|A_n|^2 = \frac{|B_n|^2}{c^2 \rho_n^2} = |G(\rho_n)|^2 \frac{\Delta \rho_n}{2\pi}. \quad (9)$$

This is the equipartition condition, with G being a deterministic variable. The equipartition condition physically means that the different modes contribute an equal amount of energy (Snieder et al. 2010). In this way, eq. (7) loses its dependence on time and becomes

$$\phi_{ZZ}(\xi) = \sum_{n=-\infty}^{\infty} |G(\rho_n)|^2 \frac{\Delta \rho_n}{2\pi} \exp(i\rho_n \xi). \quad (10)$$

As shown by Aki (1957), we can transition from the discrete wavenumber $\rho_n$ to a continuous variable $\rho$ and change the summation to an integral

$$\phi_{ZZ}(\xi) = \frac{1}{2\pi} \int_{-\infty}^{\infty} |G(\rho)|^2 \exp(i\rho\xi) d\rho. \quad (11)$$

The crux of Aki's method is that we can map the wavenumber $\rho$ integral to a frequency $\omega$ integral by a change of variables $\rho = \omega/c$

$$\phi_{ZZ}(\xi) = \frac{1}{2\pi} \int_{-\infty}^{\infty} |G(\omega/c)|^2 \exp\left(\frac{i\omega}{c}\xi\right) \frac{d\rho}{d\omega} d\omega. \quad (12)$$

Aki (1957) further argues that the term $|G(\omega/c)|^2 d\rho/d\omega$ is equal to the spectral power $\Phi(\omega)$. As a result

$$\phi_{ZZ}(\xi) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \Phi(\omega) \exp\left(\frac{i\omega}{c}\xi\right) d\omega. \quad (13)$$
For a narrowband signal centred at $\omega = \pm \omega_0$, $\Phi(\omega) = \pi P(\omega_0)\delta(\omega - \omega_0) + \pi P(-\omega_0)\delta(\omega + \omega_0)$. If $P(\omega_0) = P(-\omega_0)$, the correlation coefficient can then finally be written as

$$ \phi_{ZZ}(\xi, \omega_0) = P(\omega_0) \cos \left[ \frac{\omega_0}{c(\omega_0)} \xi \right]. \quad (14) $$

This equation shows that in 1-D the correlation coefficient is given by the cosine function (Aki 1957). This result also holds for the correlation coefficient of the horizontal components $\phi_{XX}$, except that $\phi_{XX}$ is multiplied by a factor equal to the squared ellipticity, $R^2$.

We now analyse the correlation coefficient that results from mixing vertical and horizontal components. These were not considered by Aki (1957) in his original paper. First we look at the ZX correlation

$$ \phi_{ZX}(\xi, t) = u_z(x, t) u_x(x + \xi, t). \quad (15) $$

Similar to eq. (7), we write the expression for the correlation coefficient by exploiting the fact that the different modes are uncorrelated random variables

$$ \phi_{ZX}(\xi, t) = R \sum_{n=-\infty}^{\infty} \left[ -\frac{A_n B_{-n}}{c \rho_n} \operatorname{sgn}(\rho_{-n}) \cos^2(c \rho_n t) + \frac{A_{-n} B_n}{c \rho_n} \operatorname{sgn}(\rho_{-n}) \sin^2(c \rho_n t) \right] \exp(i \rho_n \xi), \quad (16) $$

where now the equipartition condition in eq. (9) has removed terms with dependence on time as $\sin(c\rho_n t)\cos(c\rho_n t)$. We previously discussed that the condition $A_n B_{-n} = -A_{-n} B_n$ removed the terms with dependence on time as $\sin(c\rho_n t)\cos(c\rho_n t)$ from the expression for $\phi_{ZZ}$. Here, we use the relation $A_n B_{-n} = -A_{-n} B_n$ to remove the time-dependence from $\phi_{ZX}$. We arrive at

$$ \phi_{ZX}(\xi) = R \sum_{n=-\infty}^{\infty} -\frac{A_n B_{-n}}{c \rho_n} \operatorname{sgn}(\rho_{-n}) \exp(i \rho_n \xi) = R \sum_{n=-\infty}^{\infty} \frac{A_n B_{-n}}{c \rho_n} \operatorname{sgn}(\rho_n) \exp(i \rho_n \xi). \quad (17) $$

For the term containing the product of $A_n B_{-n}$ in the above equation, we substitute the following relationship

$$ \frac{A_n B_{-n}}{c \rho_n} = -\frac{A_{-n} B_n}{c \rho_n} = i |G(\rho_n)|^2 \frac{\Delta \rho_n}{2\pi}. \quad (18) $$

This relationship agrees with the equipartition condition, eq. (9), in that

$$ \frac{A_n B_{-n}}{c \rho_n} \times \frac{A_{-n} B_n}{c \rho_n} = |G(\rho_n)|^4 \left( \frac{\Delta \rho_n}{2\pi} \right)^2 \quad (19) $$

which is the same as

$$ |A_n|^2 \times \frac{|B_n|^2}{c^2 \rho_n^2} = |G(\rho_n)|^4 \left( \frac{\Delta \rho_n}{2\pi} \right)^2 \quad (20) $$

as it should be. Furthermore, the condition in eq. (18) reflects the fact that the terms $A_n B_{-n}$ and $A_{-n} B_n$ are purely imaginary. Finally, the sign convention we use in eq. (18) is consistent with the convention used to define Rayleigh wave ellipticity in Aki & Richards (1980, 2002), in which retrograde motion corresponds to a negative value for R. We further demonstrate this consistency in the following section.

Substituting eq. (18) into the summation in eq. (17) results in a Riemann sum in terms of the discrete wavenumber

$$ \phi_{ZX}(\xi) = R \sum_{n=-\infty}^{\infty} i |G(\rho_n)|^2 \frac{\Delta \rho_n}{2\pi} \operatorname{sgn}(\rho_n) \exp(i \rho_n \xi), \quad (21) $$

which, as discussed before, can be approximated in the limit as an integral

$$ \phi_{ZX}(\xi) = \frac{R}{2\pi} \int_{-\infty}^{\infty} i |G(\rho)|^2 \operatorname{sgn}(\rho) \exp(i\rho\xi) d\rho. \quad (22) $$

Changing the variable of integration to frequency instead of wavenumber gives

$$ \phi_{ZX}(\xi) = \frac{R}{2\pi} \int_{-\infty}^{\infty} i |G(\omega/c)|^2 \operatorname{sgn}(\omega/c) \exp\left(\frac{i\omega}{c}\xi\right) \frac{d\rho}{d\omega} d\omega. \quad (23) $$

Making use of the result by Aki (1957) that the term $|G(\omega/c)|^2 d\rho/d\omega$ is equal to the spectral power $\Phi(\omega)$ yields

$$ \phi_{ZX}(\xi) = \frac{R}{2\pi} \int_{-\infty}^{\infty} \Phi(\omega) \operatorname{sgn}(\omega) \exp\left(\frac{i\omega}{c}\xi\right) d\omega, \quad (24) $$

where we set $\operatorname{sgn}(\omega) = \operatorname{sgn}(\omega/c)$ since the phase velocity $c$ is always positive. For a narrowband signal centred at $\omega = \pm \omega_0$, $\Phi(\omega) = \pi P(\omega_0)\delta(\omega - \omega_0) + \pi P(-\omega_0)\delta(\omega + \omega_0)$ and, when $P(\omega_0) = P(-\omega_0)$, the correlation coefficient becomes

$$ \phi_{ZX}(\xi, \omega_0) = -P(\omega_0)R \operatorname{sgn}(\omega_0) \sin \left[ \frac{\omega_0}{c(\omega_0)} \xi \right] = -P(\omega_0)R \sin \left[ \frac{|\omega_0|}{c(\omega_0)} \xi \right]. \quad (25) $$

This equation shows that in 1-D the correlation coefficient for the ZX correlation is given by the sine function. An important aspect of this derivation is the presence of the absolute value for the frequency $|\omega_0|$ in eq. (25). As we will see, the absolute value ensures that the correlation function $\phi_{ZX}$ is closely related to the ZX component of the Rayleigh wave Green's function in 1-D. Since the imaginary parts of the two terms with $A_n$ and $B_n$ in eq. (18) are opposite in sign, the XZ component of the correlation coefficient matrix $\phi_{XZ}$ is equal to the negative of $\phi_{ZX}$.

Eq. (25) foreshadows the 2-D case since we know from Aki (1957) that the ZZ correlation is the zero order Bessel function $J_0$. By analogy between cosine and $J_0$, we can expect the ZR correlation to be dependent on the first order Bessel function $J_1$, since $J_1$ is analogous
with sine. In a later section, we will rigorously show that this is indeed the case. We can summarize our results for 1-D Rayleigh waves as follows:

$$ \phi^R(x, \omega) = \begin{bmatrix} \phi_{ZZ}^R & \phi_{ZX}^R \\ \phi_{XZ}^R & \phi_{XX}^R \end{bmatrix} = P(\omega) \begin{bmatrix} \cos(\omega x/c) & -R \sin(|\omega|x/c) \\ R \sin(|\omega|x/c) & R^2 \cos(\omega x/c) \end{bmatrix}. \quad (26) $$

In the next section, we discuss the relationship between this matrix and the Green's function matrix for Rayleigh waves in 1-D.

### 3 RELATION TO THE GREEN'S FUNCTION IN 1-D

As in Nakahara (2006), we now show that there is a relationship between the correlation matrix and the Green's function matrix. The Green's matrix for Rayleigh waves is given by eqs (7.108)–(7.110) in Aki & Richards (1980) or eqs (7.109)–(7.111) in Aki & Richards (2002). The matrix for the fundamental mode Rayleigh wave is

$$ G(x, \omega) = \begin{bmatrix} G_{ZZ} & G_{ZX} \\ G_{XZ} & G_{XX} \end{bmatrix} = \frac{1}{2cUI_1} \times \frac{c}{2\omega} \times \begin{Bmatrix} r_2^2 \exp\left(\frac{i\omega}{c}x + i\frac{\pi}{2}\right) & r_1r_2 \exp\left[\frac{i\omega}{c}x + i\frac{\pi}{2} + \operatorname{isgn}(\omega)\frac{\pi}{2}\right] \\ r_1r_2 \exp\left[\frac{i\omega}{c}x + i\frac{\pi}{2} - \operatorname{isgn}(\omega)\frac{\pi}{2}\right] & r_1^2 \exp\left(\frac{i\omega}{c}x + i\frac{\pi}{2}\right) \end{Bmatrix}, \quad (27) $$

where $r_1$ and $r_2$ are values of the horizontal and vertical components of the Rayleigh wave eigenfunction at the surface, $U$ is the group velocity and $I_1$ is an integral over depth of density weighted by the eigenfunctions at depth. Note the presence of the Hilbert transform operator $\exp[-\operatorname{isgn}(\omega)\frac{\pi}{2}] = -\operatorname{isgn}(\omega)$ in the off-diagonal components. For our purposes, we take the normalization as $r_2^2/(2cUI_1) = 1$, which is slightly different than the normalization used by Snieder (2002), $1/(8cUI_1) = 1$. With this normalization, $r_1/r_2$ is equal to the ratio of the horizontal-to-vertical motion, $R$. Thus the Green's matrix becomes

$$ G(x, \omega) = \frac{c}{2\omega} \begin{Bmatrix} \exp\left(\frac{i\omega}{c}x + i\frac{\pi}{2}\right) & R \exp\left[\frac{i\omega}{c}x + i\frac{\pi}{2} + \operatorname{isgn}(\omega)\frac{\pi}{2}\right] \\ R \exp\left[\frac{i\omega}{c}x + i\frac{\pi}{2} - \operatorname{isgn}(\omega)\frac{\pi}{2}\right] & R^2 \exp\left(\frac{i\omega}{c}x + i\frac{\pi}{2}\right) \end{Bmatrix}. \quad (28) $$

Comparing eq. (26) to eq. (28), Nakahara (2006) found the relationship between the correlation coefficient and the Green's function in 1-D to be

$$ \phi_{ZZ} = -\frac{i\omega}{c}(G_{ZZ} - G_{ZZ}^*) = \frac{2\omega}{c}\operatorname{Im}(G_{ZZ}) \quad (29) $$

for the unpolarized ZZ component. Note that, for this relation to hold, the normalization $P(\omega) = 1$ is used in eq. (26). This relationship shows that the correlation coefficient is proportional to the time-derivative of the difference between the Green's function and the time-reversed Green's function. Since $\phi_{ZZ}$ and $\phi_{XX}$ are proportional to each other, this relation holds for the XX component.

We now show that the relationship holds for the off-diagonal components as well. First, consider the ZX component for positive frequencies $\omega > 0$. From eq. (28), we find that

$$ \operatorname{Im}(G_{ZX}) = -\frac{Rc}{2\omega} \sin(\omega x/c) = -\frac{Rc}{2\omega} \sin(|\omega|x/c). \quad (30) $$

For negative frequencies $\omega < 0$, we also find

$$ \operatorname{Im}(G_{ZX}) = \frac{Rc}{2\omega} \sin(\omega x/c) = -\frac{Rc}{2\omega} \sin(|\omega|x/c). \quad (31) $$

such that, in both cases, $\phi_{ZX} = (2\omega/c) \times \operatorname{Im}(G_{ZX}) = -R\sin(|\omega|x/c)$ and the expression for $\phi_{ZX}$ in eq. (26) is verified. The expression for $\phi_{XZ}$ can similarly be verified and the connection between any element of the matrix of correlation coefficients $\phi_{ij}$ and the Green's function matrix $G_{ij}$ can be written as

$$ \phi_{ij} = -\frac{i\omega}{c}(G_{ij} - G_{ij}^*) = \frac{2\omega}{c}\operatorname{Im}(G_{ij}). \quad (32) $$

We now turn to expressions for the correlation coefficient matrix in the time domain. As shown by Nakahara (2006), the expression for $\phi_{ZZ}$ in the time domain follows as:

$$ \phi_{ZZ}(x, t) = F^{-1}[\cos(\omega x/c)] = \frac{1}{2}\delta\left(t - \frac{x}{c}\right) + \frac{1}{2}\delta\left(t + \frac{x}{c}\right), \quad (33) $$

where $F^{-1}$ denotes the inverse Fourier transform. In the context of Rayleigh waves, this result is applicable to a homogeneous half-space since $c$ generally varies as a function of frequency. However, the result is applicable to Rayleigh waves within a narrow frequency band defined by the power spectrum $P(\omega_0)$, within which $c$ does not vary greatly. For the $\phi_{ZX}$ component, we need to compute

$$ \phi_{ZX}(x, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} -R \sin(|\omega|x/c) \exp(-i\omega t) d\omega = \frac{R}{4\pi} \int_{-\infty}^{\infty} i \left[ \exp\left(\frac{i|\omega|}{c}x\right) - \exp\left(-\frac{i|\omega|}{c}x\right) \right] \exp(-i\omega t) d\omega, \quad (34) $$

where similarly to $\phi_{ZZ}$, we have assumed $R$ and $c$ are frequency independent. We proceed by breaking the integral in eq. (34) into two integrals over positive and negative frequencies

$$ \phi_{ZX}(x, t) = \frac{R}{4\pi} \int_{0}^{\infty} i \left[ \exp\left(\frac{i\omega}{c}x\right) - \exp\left(-\frac{i\omega}{c}x\right) \right] \exp(-i\omega t) d\omega + \frac{R}{4\pi} \int_{-\infty}^{0} i \left[ \exp\left(-\frac{i\omega}{c}x\right) - \exp\left(\frac{i\omega}{c}x\right) \right] \exp(-i\omega t) d\omega. \quad (35) $$
We then recombine the integrals into two integrals over the entire frequency range

$$ \phi_{ZX}(x, t) = \frac{R}{4\pi} \int_{-\infty}^{\infty} i \operatorname{sgn}(\omega) \exp\left(\frac{i\omega}{c}x\right) \exp(-i\omega t) \, d\omega - \frac{R}{4\pi} \int_{-\infty}^{\infty} i \operatorname{sgn}(\omega) \exp\left(-\frac{i\omega}{c}x\right) \exp(-i\omega t) \, d\omega. \quad (36) $$

From these expressions, it is clear that we are taking the inverse Fourier transform of the Hilbert transform operator $-i\operatorname{sgn}(\omega)$. This transform is known as the Hilbert transform $H$ of the delta function, $H[\delta(t)] = -1/\pi t$. We refer to this function in the shorthand notation $H[\delta(t)] = \delta^H(t)$. As a result, the $\phi_{ZX}$ component in the time domain can be written as

$$ \phi_{ZX}(x, t) = \frac{R}{2} \delta^H \left( t + \frac{x}{c} \right) - \frac{R}{2} \delta^H \left( t - \frac{x}{c} \right). \quad (37) $$

Whereas eq. (33) consisted of two propagating delta functions, eq. (37) indicates that there are two propagating waveforms which are the Hilbert transform of the delta function. Although the Hilbert transform of the delta function is an asymmetric wavelet, the two waveforms are subtracted and so the ZX correlation is symmetric with respect to $t = 0$. Thus, it is possible to average the positive and negative lags of the correlation in the time domain to obtain the so-called symmetric component (Lin et al. 2008). The Hilbert transform of the delta function is the waveform since, relative to $\phi_{ZZ}$, the ZX component is phase shifted by 90°. This is a direct result of the elliptical polarization of the Rayleigh wave, as exploited in Van Wijk et al. (2011).

The final issue we discuss for the 1-D case concerns the correlation coefficients when the incident energy is anisotropic. Although we did not point it out at the time, our assumption in the discussion after eq. (13), that the power spectrum was an even function of frequency (or, by $\rho = \omega/c$, an even function of wavenumber), is equivalent to assuming that the incident energy is equal from both directions. When this is not the case, we find the general form for $\phi_{ZZ}(x, \omega_0)$ to be

$$ \phi_{ZZ}(x, \omega_0) = \frac{P(\omega_0)}{2} \exp\left(\frac{i\omega_0}{c}x\right) + \frac{P(-\omega_0)}{2} \exp\left(-\frac{i\omega_0}{c}x\right), \quad (38) $$

where $P(\omega_0) \neq P(-\omega_0)$. If we set $P(\omega_0) = 2A$ and $P(-\omega_0) = 2B$ and transform to the time domain, we find, in agreement with Nakahara (2006), that

$$ \phi_{ZZ}(x, t) = A\delta\left(t - \frac{x}{c}\right) + B\delta\left(t + \frac{x}{c}\right). \quad (39) $$

In other words, the delta functions at negative and positive time-lags do not have the same amplitude. The presence of anisotropic incident energy in 1-D does not cause there to be acausal 'ghost' waveforms, as exist for anisotropic incident energy in 2-D (Nakahara 2006). By analogy with $\phi_{ZZ}$, we find that the general form for $\phi_{ZX}(x, \omega_0)$ in 1-D is given by

$$ \phi_{ZX}(x, \omega_0) = iR \left[ \frac{P(\omega_0)}{2} \operatorname{sgn}(\omega_0) \exp\left(\frac{i\omega_0}{c}x\right) + \frac{P(-\omega_0)}{2} \operatorname{sgn}(-\omega_0) \exp\left(-\frac{i\omega_0}{c}x\right) \right]. \quad (40) $$

Again setting $P(\omega_0) = 2A$ and $P(-\omega_0) = 2B$ and transforming to the time domain, we find the general form to be

$$ \phi_{ZX}(x, t) = R \left[ B \delta^H \left( t + \frac{x}{c} \right) - A \delta^H \left( t - \frac{x}{c} \right) \right]. \quad (41) $$

Just as seen with $\phi_{ZZ}$, the right- and left-going waveforms in this case are not equal in amplitude. Before leaving the 1-D case, note that the difference between the expressions in eqs (13) and (24) is the presence of the term $\operatorname{Risgn}(\omega)$ within the integrand of eq. (24). This means that $\phi_{ZX}$ is related to $\phi_{ZZ}$ by multiplication by $R$ and a Hilbert transform. We utilize this property in the 2-D case to simplify the derivation for the vertical-radial correlation $\phi_{ZR}$.

## 4 THE 2-D CASE FOR RAYLEIGH WAVES

We turn our attention to Rayleigh waves along a 2-D surface. The previous development of the 1-D case serves as a useful guide for the 2-D case. This was also the approach taken by Aki (1957) in the original derivation of SPAC. By analogy with eqs (3) and (4), the three-component displacements associated with Rayleigh waves in 2-D satisfy

$$ u_z(x, y, t) = \sum_n \sum_m A_{nm} \exp[i\rho_n r \cos(\theta_m - \psi)] \cos(c\rho_n t) + \sum_n \sum_m \frac{B_{nm}}{c\rho_n} \exp[i\rho_n r \cos(\theta_m - \psi)] \sin(c\rho_n t) \quad (42) $$

$$ u_r(x, y, t) = R \sum_n \sum_m -A_{nm} \cos(\theta_m - \psi) \exp[i\rho_n r \cos(\theta_m - \psi)] \sin(c\rho_n t) \\ + R \sum_n \sum_m \frac{B_{nm}}{c\rho_n} \cos(\theta_m - \psi) \exp[i\rho_n r \cos(\theta_m - \psi)] \cos(c\rho_n t) \quad (43) $$

$$ u_{\psi}(x, y, t) = R \sum_{n} \sum_{m} -A_{nm} \sin(\theta_m - \psi) \exp[i\rho_n r \cos(\theta_m - \psi)] \sin(c\rho_n t) \\ + R \sum_{n} \sum_{m} \frac{B_{nm}}{c\rho_n} \sin(\theta_m - \psi) \exp[i\rho_n r \cos(\theta_m - \psi)] \cos(c\rho_n t), \quad (44) $$

where $r = \sqrt{x^2 + y^2}$ and $\psi = \tan^{-1}(y/x)$ are the radial distance and azimuthal angle describing the relative positions of the two seismometers. In contrast to the 1-D case, the displacements are now the result of two summations over the radius $\rho_n$ and polar angle $\theta_m$ for the wavenumber.
Another difference with the 1-D case is that the factor sgn($\rho_n$) does not appear in eqs (43) and (44) as it did in eq. (4). This is because the sign change due to different propagation directions is implicitly handled by the factors cos($\theta_m$ - $\psi$) and sin($\theta_m$ - $\psi$). Finally, note that the Rayleigh waves in eqs (42)-(44) do not necessarily propagate along the direct path between the receivers and, as a result, the horizontal particle motion is in general not strictly radial.

We first consider the spatial correlation of the vertical components (ZZ) before moving on to the mixed-component ZR correlation

$$ \phi_{ZZ}(r, \psi, t) = \overline{u_z(x_0, y_0, t) u_z(x_0 + r \cos \psi, y_0 + r \sin \psi, t)}. \quad (45) $$

Similar to the 1-D case, we find that

$$ \phi_{ZZ}(r, \psi, t) = \sum_n \sum_m \left[ |A_{nm}|^2 \cos^2(c\rho_n t) + \frac{|B_{nm}|^2}{c^2 \rho_n^2} \sin^2(c\rho_n t) \right] \exp[i\rho_n r \cos(\theta_m - \psi)]. \quad (46) $$

The equipartition condition in 2-D relates the $A_{nm}$ and $B_{nm}$ coefficients as

$$ |A_{nm}|^2 = \frac{|B_{nm}|^2}{c^2 \rho_n^2} = |G(\rho_n, \theta_m)|^2 \frac{\rho_n \Delta\rho_n \Delta\theta_m}{(2\pi)^2}. \quad (47) $$

Inserting this into eq. (46), we obtain the following expression

$$ \phi_{ZZ}(r, \psi) = \sum_n \sum_m |G(\rho_n, \theta_m)|^2 \frac{\rho_n \Delta\rho_n \Delta\theta_m}{(2\pi)^2} \exp[i\rho_n r \cos(\theta_m - \psi)]. \quad (48) $$

As shown previously in 1-D, we can transition from the discrete wavenumber $\rho_n$ and discrete polar angle $\theta_m$ to the continuous variables $\rho$ and $\theta$. This transition changes the summation to an integral

$$ \phi_{ZZ} = \frac{1}{(2\pi)^2} \int_0^{\infty} \int_0^{2\pi} |G(\rho, \theta)|^2 \exp[i\rho r \cos(\theta - \psi)] \rho \, d\rho \, d\theta. \quad (49) $$

From Cox (1973), the exponential in this equation can be expressed as a summation over cylindrical harmonics

$$ \exp[i\rho r \cos(\theta - \psi)] = \sum_{m=-\infty}^{\infty} i^m J_m(\rho r) \exp[-im(\theta - \psi)]. \quad (50) $$

Inserting this infinite series into eq. (49) gives

$$ \phi_{ZZ} = \frac{1}{2\pi} \int_{0}^{\infty} \sum_{m=-\infty}^{\infty} i^m J_m(\rho r) \left[ \frac{1}{2\pi} \int_{0}^{2\pi} |G(\rho, \theta)|^2 \exp[-im(\theta - \psi)] d\theta \right] \rho \, d\rho. \quad (51) $$

We further assume that $G$ can be split into separate terms dependent on $\rho$ and $\theta$ as $|G(\rho, \theta)|^2 = |G(\rho)|^2 p(\theta)$. Inserting this into the above expression results in

$$ \phi_{ZZ} = \frac{1}{2\pi} \int_{0}^{\infty} |G(\rho)|^2 \rho \sum_{m=-\infty}^{\infty} i^m J_m(\rho r) \gamma_m \, d\rho, \quad (52) $$

where the angular variation determines the following factor $\gamma_m$

$$ \gamma_m = \frac{1}{2\pi} \int_{0}^{2\pi} p(\theta) \exp[-im(\theta - \psi)] d\theta. \quad (53) $$

As in the 1-D case, we implement the change of variables $\rho = \omega/c$ to obtain

$$ \phi_{ZZ} = \frac{1}{\pi} \int_{0}^{\infty} \frac{1}{2} |G(\omega/c)|^2 \frac{\omega}{c} \frac{d\rho}{d\omega} \sum_{m=-\infty}^{\infty} i^m J_m \left( \frac{\omega}{c} r \right) \gamma_m \, d\omega. \quad (54) $$

As shown by Aki (1957), we find that the spectral power is given by

$$ \Phi(\omega) = \frac{1}{2} |G(\omega/c)|^2 \frac{\omega}{c} \frac{d\rho}{d\omega} \, d\omega \quad (55) $$

and upon substituting this expression into eq. (54), we find that

$$ \phi_{ZZ} = \frac{1}{\pi} \int_{0}^{\infty} \Phi(\omega) \sum_{m=-\infty}^{\infty} i^m J_m \left( \frac{\omega}{c} r \right) \gamma_m \, d\omega. \quad (56) $$

Finally, we take the wavefield to be narrowband by assuming $\Phi(\omega) = \pi P(\omega_0) \delta(\omega - \omega_0)$, which yields

$$ \phi_{ZZ} = P(\omega_0) \sum_{m=-\infty}^{\infty} i^m J_m \left( \frac{\omega_0}{c} r \right) \gamma_m. \quad (57) $$

Since $p(\theta)$ is real, $\gamma_{-m} = \gamma_m^*$ and the two-sided infinite series can be rewritten as

$$ \phi_{ZZ} = P(\omega_0) \sum_{m=0}^{\infty} \epsilon_m i^m J_m \left( \frac{\omega_0}{c} r \right) \text{Re}[\gamma_m] \quad (58) $$

where $\epsilon_0 = 1$ and $\epsilon_m = 2$ for $m \neq 0$. For an isotropic wavefield, one in which energy propagates equally in all directions, $p(\theta) = 1$, $\gamma_0 = 1$ and $\gamma_m = 0$ for $m \neq 0$. In this case,
$$ \phi_{ZZ} = P(\omega_0)J_0\left(\frac{\omega_0}{c}r\right) \qquad (59) $$

which is the well-known result of Aki (1957) for unpolarized waves.

Using the methodology demonstrated above for the ZZ correlation, we now consider the ZR correlation given by

$$ \phi_{ZR}(r, \psi, t) = \frac{u_z(x_0, y_0, t) u_r(x_0 + r \cos \psi, y_0 + r \sin \psi, t)}{u_z(x_0, y_0, t) u_z(x_0 + r \cos \psi, y_0 + r \sin \psi, t)}. \qquad (60) $$

The derivation in this case follows similar logic as the above derivation for $\phi_{ZZ}$ and therefore we only point out the differences for $\phi_{ZR}$. Due to the term $\cos(\theta - \psi)$ in eq. (43), the factor $\gamma_m$ is different from the $\phi_{ZZ}$ case and is given by

$$ \gamma_m = \frac{1}{2\pi} \int_0^{2\pi} p(\theta) \cos(\theta - \psi) \exp[-i m(\theta - \psi)] d\theta. \qquad (61) $$

As pointed out in the 1-D case, $\phi_{ZR}$ is different from $\phi_{ZZ}$ in that the term $R$ isgn($\omega$) acts as an additional frequency filter for $\phi_{ZR}$. This results in the following expression for $\phi_{ZR}$

$$ \phi_{ZR} = \frac{R}{\pi} \int_{0}^{\infty} i \operatorname{sgn}(\omega) \Phi(\omega) \sum_{m=-\infty}^{\infty} i^m J_m \left(\frac{\omega}{c}r\right) \gamma_m \, d\omega. \qquad (62) $$

Again, we take $\Phi(\omega) = \pi P(\omega_0) \delta(\omega - \omega_0)$, resulting in

$$ \phi_{ZR} = R P(\omega_0) i \operatorname{sgn}(\omega_0) \sum_{m=-\infty}^{\infty} i^m J_m \left(\frac{\omega_0}{c} r \right) \gamma_m. \qquad (63) $$

Since $p(\theta)$ is real, $\gamma_{-m} = \gamma_m^*$ and the two-sided infinite series can be rewritten as

$$ \phi_{ZR} = R P(\omega_0) i \operatorname{sgn}(\omega_0) \sum_{m=0}^{\infty} \epsilon_m i^m J_m \left(\frac{\omega_0}{c} r \right) \operatorname{Re}[\gamma_m], \qquad (64) $$

where as before, $\epsilon_0 = 1$ and $\epsilon_m = 2$ for $m \neq 0$.

For an isotropic wavefield, $p(\theta) = 1$, $\gamma_1 = 1/2$ and $\gamma_m = 0$ for $m \neq 1$. In this case,

$$ \phi_{ZR} = -R P(\omega_0) J_1\left(\frac{\omega_0}{c} r\right) \operatorname{sgn}(\omega_0) = -R P(\omega_0) J_1\left(\frac{|\omega_0|}{c} r\right), \qquad (65) $$

where the absolute value follows from the fact that $J_m(x)$ is an odd function of $x$ if $m$ is an odd number, with $J_m(x)$ an even function otherwise. This result shows that the correlation function for the ZR component is given by the first-order Bessel function $J_1$ in the case of an isotropic wavefield as opposed to the zeroth-order Bessel function for the ZZ component. Although results for SPAC are often given as normalized correlation functions, with the normalization at $r = 0$ (Aki 1957; Chouet 1996), here we can see that such a normalization is not possible for the ZR correlation since $\phi_{ZR}$ is equal to 0 at $r = 0$. As a result, we have chosen to present correlation functions that are not normalized throughout this paper. The lack of normalization means that all of our expressions contain the spectral power at the reference frequency $P(\omega_0)$.

As we saw previously for the 1-D case, it can be shown that the RZ correlation is the negative of the ZR correlation. Although the two correlations are equal to each other in magnitude, they do constitute independent measurable quantities since they utilize different components of the wavefield. Thus it is correct to consider both the ZR and RZ correlations as independent quantities. The RR and TT correlations were considered previously by Aki (1957) for an isotropic Rayleigh wavefield and shown to be given by

$$ \phi_{RR} = R^2 P(\omega) [J_0(\omega r/c_R) - J_2(\omega r/c_R)]/2, \qquad (66) $$

and

$$ \phi_{TT} = R^2 P(\omega) [J_0(\omega r/c_R) + J_2(\omega r/c_R)]/2. \qquad (67) $$

Note that Aki (1957) used normalized correlation functions whereas these are not normalized. A surprising result from Aki (1957) is that the TT correlation in the presence of an isotropic Rayleigh wavefield is not zero. It is important to remember that asymptotically $J_2(x) = -J_0(x)$ for large $x$. Thus, in the far field, $\phi_{TT} = 0$ and $\phi_{RR}$ is proportional to $\phi_{ZZ}$ by a factor of $R^2$, that is $\phi_{RR} = R^2 P(\omega) J_0(\omega r/c_R)$. However, in the near field, $\phi_{TT}$ is non-zero. Due to the connections with the time-domain Green's function, as shown in the next section, this means that in the near-field of a transverse force source, Rayleigh waves are measured on a transverse receiver.

The remaining correlations involve combinations of Z or R components with the transverse component: $ZT$, $TZ$, $RT$ and $TR$. For an isotropic wavefield, all of these terms are equal to zero. In Appendix A, we show the $\gamma_m$ for these components in the case of generally anisotropic incidence. For isotropic incidence, the $\gamma_m$ for these components are given in Appendix B by setting the angular range variable $\Delta = 180^{\circ}$. Following this procedure for the $ZT$, $TZ$, $RT$ and $TR$ components shows that they are indeed equal to zero for isotropic incidence.

As shown previously, we can summarize the results for Rayleigh waves in 2-D in the presence of isotropic incidence with the complete correlation matrix

$$ \phi^R(r, \omega) = \begin{bmatrix} \phi_{ZZ}^R & \phi_{ZR}^R & \phi_{TZ}^R \\ \phi_{RZ}^R & \phi_{RR}^R & \phi_{RT}^R \\ \phi_{TZ}^R & \phi_{TR}^R & \phi_{TT}^R \end{bmatrix} = P^R(\omega) \times \begin{bmatrix} J_0(\omega r/c_R) & -R J_1(|\omega r/c_R|) & 0 \\ R J_1(|\omega r/c_R|) & R^2 [J_0(\omega r/c_R) - J_2(\omega r/c_R)]/2 & 0 \\ 0 & 0 & R^2 [J_0(\omega r/c_R) + J_2(\omega r/c_R)]/2 \end{bmatrix}, \qquad (68) $$
where $P^R(\omega)$ is the power spectrum of the Rayleigh waves. As shown in Appendix A, anisotropic incidence generally causes the off diagonal terms that are zero for isotropic incidence to become non-zero. Anisotropic incidence also causes the correlation matrix to become complex-valued. For Love waves, all off-diagonal terms are zero in the case of isotropic incidence, meaning that only RR and TT components exist since ZZ is equal to zero. Thus, the results of Aki (1957) for the main diagonal are complete for Love waves. Note that, in the presence of both Rayleigh and Love wave energy, the ZZ, ZR and RZ components are not affected by the Love waves. This is in contrast to the RR and TT coefficients, which are sensitive to both Rayleigh and Love waves, at least within the near-field as described above. In Appendix A, we give expressions for both Rayleigh and Love wave correlation matrices in the case of generally anisotropic incidence.

## 5 RELATION TO THE GREEN'S FUNCTION IN 2-D

As in the 1-D case, we now show that there is a relationship between the 2-D correlation matrix and the Green's function matrix. The Green's matrix for Rayleigh waves is given in the far-field by eq. (7.146) in Aki & Richards (1980) or eq. (7.147) in Aki & Richards (2002). The matrix for the fundamental mode Rayleigh wave is

$$ G(x, \omega) = \begin{bmatrix} G_{ZZ} & G_{ZR} \\ G_{RZ} & G_{RR} \end{bmatrix} = \frac{1}{8cU I_1} \times \sqrt{\frac{2c}{\pi |\omega|r}} \times \left\{ \begin{array}{ll} r_2^2 \exp\left[\frac{i\omega}{c}r + i \operatorname{sgn}(\omega)\frac{\pi}{4}\right] & r_1 r_2 \exp\left[\frac{i\omega}{c}r + i \operatorname{sgn}(\omega)\frac{3\pi}{4}\right] \\ r_1 r_2 \exp\left(\frac{i\omega}{c}r - i \operatorname{sgn}(\omega)\frac{\pi}{4}\right) & r_1^2 \exp\left[\frac{i\omega}{c}r + i \operatorname{sgn}(\omega)\frac{\pi}{4}\right] \end{array} \right\}, \quad (69) $$

where $r_1$ and $r_2$ are values of the radial and vertical components of the Rayleigh wave eigenfunction at the surface, $U$ is the group velocity and $I_1$ is an integral over depth of density weighted by the eigenfunctions at depth. Note that the Green's matrix in Aki & Richards (1980, 2002) is a $3 \times 3$ in terms of the three Cartesian coordinates $x$, $y$ and $z$. Here, we use a $2 \times 2$ matrix in terms of the radial and vertical components, $r$ and $z$, since the Rayleigh waves only exist on those components in the far-field. We again take the normalization as $r_2^2/(2cU I_1) = 1$. With this normalization, $r_1/r_2$ is equal to the ratio of the horizontal-to-vertical motion, $R$. Thus the Green's matrix becomes

$$ G(r, \omega) = \frac{1}{4} \times \sqrt{\frac{2c}{\pi |\omega|r}} \left\{ \begin{array}{ll} \exp\left[\frac{i\omega}{c}r + i \operatorname{sgn}(\omega)\frac{\pi}{4}\right] & R \exp\left[\frac{i\omega}{c}r + i \operatorname{sgn}(\omega)\frac{3\pi}{4}\right] \\ R \exp\left[\frac{i\omega}{c}r - i \operatorname{sgn}(\omega)\frac{\pi}{4}\right] & R^2 \exp\left(\frac{i\omega}{c}r + i \operatorname{sgn}(\omega)\frac{\pi}{4}\right) \end{array} \right\}. \quad (70) $$

Nakahara (2006) found a relationship between the correlation coefficient and the Green's function in 2-D. For the unpolarized ZZ correlation, the relation between the two is

$$ H[\phi_{ZZ}] = -2(G_{ZZ} - G_{ZZ}^*) \quad (71) $$

where $H$ denotes the Hilbert transform operator $\exp[-\operatorname{sgn}(\omega)\frac{\pi}{2}] = -\operatorname{sgn}(\omega)$. As in the 1-D case, the normalization $P(\omega) = 1$ is used for $\phi_{ZZ}$ in this expression. We point out that there is a difference in sign between eq. (71) and the result of Nakahara (2006). This stems from the fact that there is a sign difference in the definition of the 2-D Green's function in eq. (18.23) of Snieder (2001), on which the result of Nakahara (2006) was based, and the Rayleigh wave Green's function in eq. (7.100) of Aki & Richards (1980). The sign difference arises because the forcing terms shown in these two references are on opposite sides of the 2-D wave equation. Note that there was no difference in sign for the previously considered 1-D case, eq. (29), and the corresponding result of Nakahara (2006). This is because, in addition to the sign difference, a sign error exists in the expression for the 1-D Green's function in eqs (18.55) and (18.56) of Snieder (2001).

In the following, we show that the relation in eq. (71) holds for the ZR component in 2-D. In contrast to the 1-D case, we rely on asymptotic approximations of Bessel functions to show this and therefore the demonstration applies to the far-field. This is necessary since the Green's matrix eq. (69) is derived in the far-field. From Snieder (2001), we know that the Bessel function of any order can be written as the sum of Hankel functions of the same order of both the first and second kind

$$ J_m\left(\frac{\omega}{c}r\right) = \frac{1}{2}\left[H_m^{(1)}\left(\frac{\omega}{c}r\right) + H_m^{(2)}\left(\frac{\omega}{c}r\right)\right]. \quad (72) $$

For large values of $\omega r/c$, this can be expressed in the far-field as the sum of exponentials

$$ J_m\left(\frac{\omega}{c}r\right) \sim \frac{1}{2}\sqrt{\frac{2c}{\pi|\omega|r}} \left\{ \exp\left[\frac{i\omega}{c}r - i\left(m + \frac{1}{2}\right)\operatorname{sgn}(\omega)\frac{\pi}{2}\right] + \exp\left[-\frac{i\omega}{c}r + i\left(m + \frac{1}{2}\right)\operatorname{sgn}(\omega)\frac{\pi}{2}\right] \right\}. \quad (73) $$

For $\omega > 0$, we use this expression for $m = 1$ and find that

$$ -R J_1\left(\frac{|\omega_0|}{c}r\right) = -\frac{R}{2}\sqrt{\frac{2c}{\pi|\omega|r}}\left[i\exp\left(-\frac{i\omega}{c}r\right) - \exp\left(\frac{i\omega}{c}r\right)\right]\exp\left(\frac{i\pi}{4}\right). \quad (74) $$

For the Green's function, we take the following expression from the matrix in eq. (70) for $\omega > 0$

$$ G_{ZR} = \frac{1}{4} \times \sqrt{\frac{2c}{\pi |\omega|r}} R \exp\left(\frac{i\omega}{c}r + i\frac{3\pi}{4}\right) \quad (75) $$

to find that

$$ -2(G_{ZR} - G_{ZR}^*) = -\frac{R}{2}\sqrt{\frac{2c}{\pi|\omega|r}}\left[i\exp\left(\frac{i\omega}{c}r\right) + \exp\left(-\frac{i\omega}{c}r\right)\right]\exp\left(\frac{i\pi}{4}\right). \quad (76) $$
Eqs (74) and (76) can be made equal to each other if eq. (74) is multiplied by a factor of $-i$. This factor is the Hilbert transform operator $-i\mathbf{\sigma}$ for the case we are currently considering, $\omega > 0$. The equivalence we have just shown for $\omega > 0$ is demonstrated for $\omega < 0$ in Appendix C. Therefore, we have shown that

$$H[\phi_{ij}] = -2(G_{ij} - G_{ij}^*) \qquad (77)$$

holds for all components of the correlation matrix and the Green's function matrix.

## 6 RESILIENCY TO GHOSTS

The ZR and RZ correlations have advantageous properties compared to the ZZ correlation in the presence of anisotropic ambient noise. We demonstrate these properties theoretically here, but this has also been shown observationally by Van Wijk et al. (2011). To investigate the robustness of the ZR correlation, we proceed in the manner of Nakahara (2006) and study the theoretical time-domain correlation functions within a cone of times prior to the arrival of the direct wave. Since a relationship exists between correlation functions and the Green's function, there should ideally be no high-amplitude arrivals prior to the direct wave. As we show here, such arrivals—called ‘ghost’ arrivals—do occur when the noise wavefield is highly anisotropic in the out-of-plane direction between two seismometers. The ghost arrivals are more pronounced in ZZ correlations compared to ZR correlations and as a result, the ZR correlations can be considered to be more robust.

The situation we consider is shown in Fig. 1, where two stations are subjected to highly directional seismic noise centred about the azimuth $\phi_0$ and spanning an angular range of $2\Delta$. From Nakahara (2006), the time-domain expressions for the ZZ correlation and its Hilbert transform are given for $|t| < r/c$ as

$$\phi_{ZZ}(r, t) = \frac{1}{\pi} \sum_{m=0}^{m=\infty} \epsilon_m \mathrm{Re} [\gamma_m] \frac{T_m(ct/r)}{\sqrt{(r/c)^2 - t^2}} \qquad (78)$$

and

$$H[\phi_{ZZ}(r, t)] = \frac{2c}{\pi r} \sum_{m=1}^{m=\infty} \mathrm{Re}[\gamma_m] U_{m-1}\left(\frac{ct}{r}\right), \qquad (79)$$

where we take the azimuth between the stations to be $0^\circ$. In these expressions, $T_m(ct/r)$ and $U_{m-1}(ct/r)$ are the $m$th order Chebyshev polynomials of the first kind and $(m-1)$th order Chebyshev polynomials of the second kind, respectively. Nakahara (2006) showed that for a restricted angular range $\Delta$ of incident energy in a cone centred about azimuth $\phi_0$, the coefficient $\gamma_m$ is given by

$$\mathrm{Re}[\gamma_m] = \cos[m(\phi_0 - \psi)] \sin c(\Delta m). \qquad (80)$$

This result is also given in eq. (B1) in Appendix B along with other angular integrals for incident energy within a range of angles.

For the expression for the ZR correlation in the time domain, we take advantage of the fact that the main differences with the ZZ correlation are that the ZR correlation is

(i) Hilbert transformed compared to the ZZ correlation.

(ii) Multiplied by the ratio of the radial-to-vertical motion $R$.

(iii) Subject to a different $\gamma_m$ coefficient given by eq. (61) instead of eq. (53).

Taking these differences into account, we find that the time-domain expressions for the ZR correlation and its Hilbert transform are given for $|t| < r/c$ as

$$\phi_{ZR}(r, t) = -\frac{2cR}{\pi r} \sum_{m=1}^{m=\infty} \mathrm{Re}[\gamma_m] U_{m-1}\left(\frac{ct}{r}\right) \qquad (81)$$

and

$$H[\phi_{ZR}(r, t)] = \frac{R}{\pi} \sum_{m=0}^{m=\infty} \epsilon_m \mathrm{Re}[\gamma_m] \frac{T_m(ct/r)}{\sqrt{(r/c)^2 - t^2}} \qquad (82)$$

<figure><img src="image_2.png" /><figcaption>Figure 1. Angles used to demonstrate the effect of anisotropic incidence on correlations.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison of ghosts in ZZ and ZR correlations for $\phi_0 = 75^\circ$ and four types of directional noise: isotropic incidence ($\Delta = 180^\circ$), incidence from a half-plane ($\Delta = 90^\circ$), a cone of angles given by $\Delta = 30^\circ$, and a narrow cone of $\Delta = 5^\circ$. The correlations $\phi_{ZZ}$ and $\phi_{ZR}$ are given as dashed lines, with their Hilbert transforms shown as solid lines. Time has been normalized by distance and phase velocity as in Nakahara (2006). A normalized time of unity is therefore equal to the arrival time of the direct wave between the stations.</figcaption></figure>

with the coefficient $\gamma_m$ derived from eq. (B2) in Appendix B and given as

$$ \mathrm{Re}[\gamma_m] = \frac{\cos[(m-1)(\phi_0 - \psi)] \sin[(m-1)\Delta]}{2} + \frac{\cos[(m+1)(\phi_0 - \psi)] \sin[(m+1)\Delta]}{2}. \quad (83) $$

In Fig. 2, we evaluate and plot the theoretical expressions for the ZZ and ZR correlations and their Hilbert transforms in eqs (78), (79), (81) and (82) by taking terms in the summations up to $m = 10 000$. The plots are shown for a noise azimuth of $\phi_0 = 75^\circ$ and four types of directional noise: isotropic incidence ($\Delta = 180^\circ$), incidence from a half-plane ($\Delta = 90^\circ$), a cone of angles given by $\Delta = 30^\circ$ and a narrow cone of $\Delta = 5^\circ$. Comparison of the ZZ and ZR correlations shows that the ZR correlation is less susceptible than ZZ to acausal ghost arrivals that develop for $\Delta = 30^\circ$ and $\Delta = 5^\circ$. The resiliency to ghosts for ZR correlations is due to the $\gamma_m$ coefficient in eq. (61). The polarization of the waves in the ZR case acts as a spatial filter that dampens the out-of-plane ghost arrivals. This makes the ZR correlations more desirable when considering relatively short time windows for the correlations, when the condition of isotropic incidence has a smaller likelihood of being satisfied. In practice, the length of the time window in time-lapse noise correlation studies involves a tradeoff between convergence to the ideal Green's function and the desire to have adequate time resolution. For example, in Haney (2009), correlations of infrasound noise averaged over 10 hr were able to detect day-long variations in sound speed due to changing atmospheric conditions. The robustness demonstrated in Fig. 2 means the ZR correlations are a better candidate for use in detecting time-lapse changes with seismic noise than ZZ correlations. Fig. 3 shows the case when the waves propagate perpendicular to the line of receivers, when $\phi_0 = 90^\circ$. The dampening effect of the acausal ghost arrivals is observed to be even more powerful for ZR compared to ZZ. Note that in the presence of both Rayleigh and Love wave noise, the ZZ and ZR correlations are unaffected by the Love waves, even in the near field when the two stations are within a wavelength of each other. Therefore, these results hold for general ambient noise conditions, when the noise is comprised of both Rayleigh and Love waves.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Comparison of ghosts in ZZ and ZR correlations for $\phi_0 = 90^\circ$ and 4 types of directional noise: isotropic incidence ($\Delta = 180^\circ$), incidence from a half-plane ($\Delta = 90^\circ$), a cone of angles given by $\Delta = 30^\circ$, and a narrow cone of $\Delta = 5^\circ$. The correlations $\phi_{ZZ}$ and $\phi_{ZR}$ are given as dashed lines, with their Hilbert transforms shown as solid lines. Time has been normalized by distance and phase velocity as in Nakahara (2006). A normalized time of unity is therefore equal to the arrival time of the direct wave between the stations.</figcaption></figure>

# 7 SPAC WITH ZR CORRELATIONS AT AKUTAN VOLCANO

Recently, Ekström *et al.* (2009) demonstrated the connection between the SPAC method and time-domain correlations of ambient noise between vertical components (ZZ). Namely, the zeros of the real part of the spectrum can be associated with the zeros of the zeroth-order Bessel function. Here, we show that our results allow the same methodology to be applied to ZR correlations as well, by associating the zeros of the real part of the spectrum with the zeros of the first-order Bessel function.

We use data from two broadband, three-component seismometers located at Akutan Volcano, historically one of the most active volcanoes in the Aleutian Island chain (Lu *et al.* 2000). The two sensors, AKMO and AKLV, are both Guralp CMG-6TD instruments (Dixon *et al.* 2008) that exist within the network operated by the Alaska Volcano Observatory. The location of Akutan Island and the seismic stations are given in Fig. 4. We compute time-domain correlations between all components (Z, N and E) over 40 d beginning on 2007 November 1. The correlations are done in a manner that maintains the relative amplitudes between the three-components associated with a single station (Haney *et al.* 2011). This is accomplished by applying identical time-normalization and frequency-normalization filters to each of the three components. For the time-normalization, a moving average window of rms amplitude is computed at each time sample for each of the three components associated with a single station. We use a time-window of 10 s for the moving averages. A single time-normalization filter is formed from the maximum of the three moving-averages at each time sample. The three-component seismograms are then divided by the same value of the time-normalization filter at each time sample. Note that this balancing of amplitude in time normalization is impossible if the sign-bit type of normalization is applied to each of the three components individually. Hour-long time-normalized seismograms are then subjected to frequency normalization by forming smoothed versions of the amplitude spectra for each of the three components. The smoothing is done over a frequency range of 0.025 Hz. A frequency-normalization filter is formed from the maximum of the three smoothed amplitude spectra at each discrete frequency between 0.1 and 0.7 Hz. The spectra for each of the three components are then divided by the same frequency-normalization filter at each discrete frequency. Outside of this frequency band, the frequency spectrum of the signals are
<figure><img src="image_5.png" /><figcaption>**Figure 4.** A regional map showing the location of Akutan Volcano (upper panel) and a local map of Akutan Island (lower panel). The locations of the two stations used to compute cross correlations are given on the local map. Bold contours are plotted at 0 and 850 m elevation; light contours are shown between 150 and 675 m elevation in intervals of 175 m. The highest contour line, at 850 m, clearly marks the location of the summit caldera at Akutan Volcano.</figcaption></figure>

tapered to zero. Once the signals have been time and frequency normalized, all combinations of hour-long, three-component seismograms are correlated between the stations. Once we compute the correlation matrix between the stations in the (Z, N, E) coordinate system, we transform to (Z, R, T) by applying a tensor rotation

$$ \begin{bmatrix} ZZ & ZR & ZT \\ RZ & RR & RT \\ TZ & TR & TT \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \theta & \sin \theta \\ 0 & -\sin \theta & \cos \theta \end{bmatrix} \times \begin{bmatrix} ZZ & ZE & ZN \\ EZ & EE & EN \\ NZ & NE & NN \end{bmatrix} \times \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \theta & -\sin \theta \\ 0 & \sin \theta & \cos \theta \end{bmatrix}, \qquad (84) $$

where $\theta$ is the azimuth (in degrees counter-clockwise from east) between the stations. The importance of this type of processing is that the angular rotation can be performed after the cross-correlation step. This means that for an array of stations, the rotations into the (Z, R, T) coordinate system can be done on-the-fly after computing the correlations. Also note that, had the correlations not been done in such a way as to balance the three-component amplitudes, the rotation to the (Z, R, T) coordinate system would not have been accurate.

In Fig. 5, we plot the real part of the frequency spectrum of the time-domain correlations for both ZZ and ZR. It is clear that the zeros of ZZ and ZR do not coincide. This is expected since the zeros of the ZZ spectrum are associated with the zeros of the zeroth-order Bessel function, with the zeros of the ZR spectrum given by the first-order Bessel. Over the frequency band from 0.1 to 0.7, we pick four zeros in both the ZZ and ZR spectra, as shown in Fig. 5. We compute the phase velocities c at these zeros by setting the argument of the Bessel functions equal to their kth zero

$$ \frac{\omega_k r}{c} = z_k, \qquad (85) $$

where r = 8.7 km is the interstation distance and $z_k$ is the kth zero of the zeroth or first order Bessel function, depending on whether ZZ or ZR is being analysed. It is clear from Fig. 5 that the lowermost frequency, 0.1 Hz, is low enough to capture the first zero crossing on the vertical component. For greater interstation distances, the lowermost frequency must be decreased in order to detect the first zero crossing. If seismic noise does not exist for lower frequencies, there can be some ambiguity in the association of zero-crossings with a particular zero.

We plot the computed phase velocities at Akutan in Fig. 6 along with the Rayleigh wave dispersion curve for a 1-D model of Okmok Volcano developed by Masterlark et al. (2010). The phase velocities at Akutan are in good agreement with the Okmok model, showing that the values are generally consistent with those obtained at a nearby volcano in the Aleutian Arc. Uncertainty estimates are also plotted for the phase velocities in Fig. 6. The errors are obtained by fitting lines to the ZZ and ZR spectra in the immediate vicinity of the zero crossings. Close to the zero-crossings, the spectra should be almost linear, just as a sine function is linear near its zero-crossings. From the linear fit, we obtain error estimates (variances) for the slope and the y-intercept directly as well as the covariance between the slope and y-intercept
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The real part of the frequency spectrum for the average time-domain correlation between stations AKMO and AKLV at Akutan Volcano for ZZ (top panel) and ZR (bottom panel). The spectrum is tapered to zero below 0.1 Hz. The zeros of the real part of the frequency spectrum are shown for both ZZ and ZR by vertical lines. Rayleigh wave phase velocities, measured at the locations of the zero crossings of ZZ and ZR, are shown in Fig. 6.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Rayleigh wave phase velocities at Akutan Volcano measured from the first four ZZ (red squares) and ZR (red circles) zero-crossings. The phase velocities are obtained by associating the ZZ zero-crossings with $J_0$ and the ZR zero-crossings with $J_1$. Uncertainties in the phase velocities are also plotted as described in the text. The blue curve is the Rayleigh wave phase velocity dispersion curve for a 1-D model at nearby Okmok Volcano (Masterlark et al. 2010). The Okmok dispersion curve shows that the phase velocities obtained at Akutan volcano are consistent with structure inferred at a nearby volcano in the Aleutian Arc.</figcaption></figure>

estimates. Finally, the errors for the x-intercept (i.e. zero-crossing), indicated by the symbol $p$, are estimated from the variances and covariance for the slope $m$ and y-intercept $b$ by simple propagation of errors

$$ \sigma_p^2 = \sigma_b^2 \left(\frac{\partial p}{\partial b}\right)^2 + \sigma_m^2 \left(\frac{\partial p}{\partial m}\right)^2 + 2 \operatorname{cov}_{mb} \left(\frac{\partial p}{\partial b}\right) \left(\frac{\partial p}{\partial m}\right), \quad (86) $$

where the derivatives are computed according to the relation $p = -b/m$.

In Fig. 6, the phase velocity estimations from both ZZ and ZR correlations are overall consistent with each other and show a decreasing phase velocity with increasing frequency. Note how the values of phase velocity from the ZR correlations naturally interpolate the values from the ZZ correlation. The agreement between the ZZ and ZR measurements below 0.4 Hz supports the theoretical results derived in this paper. Above 0.4 Hz, there is some discrepancy between the phase velocities estimated from the ZZ and ZR correlations as the estimates from the ZR correlations appear to be slightly larger. This is likely related to the wavelength of the Rayleigh waves at 0.4 Hz approaching the interstation distance. For a nominal phase velocity of 2.5 km s⁻¹, the wavelength at 0.4 Hz is 6.25 km which is comparable to the interstation distance of 8.7 km. The assumption that the subsurface is an effective homogeneous medium between the stations becomes progressively less applicable as the frequency increases and the wavelength becomes smaller than the interstation distance. At higher frequencies, ray-bending
and scattering may also become more important. Ray-bending, in particular, poses a problem since the ZR correlation assumes the azimuth and backazimuth of Rayleigh wave propagation between the stations are directed along the interstation azimuth. Another possible cause for the discrepancy at frequencies above 0.4 Hz is if the source distribution is not close to being isotropic and the zeros of the ZZ and ZR correlations are not determined simply from the first order term, either $J_0$ for ZZ or $J_1$ for ZR. These issues deserve further attention in future work.

Additional challenges for future work would involve fitting the observed spectra over the entire frequency band for a source power distribution and a generally anisotropic incident wavefield instead of only using the zeros of the spectra. Alternatively, the $J_1$ dependence would be evident for the ZR correlation within a small-aperture array of several seismometers in which the spatial sampling was dense enough to observe the correlation coefficient as a function of interstation distance.

## 8 CONCLUSION

We have extended the SPAC method to the full matrix of correlation coefficients that exists between three-component seismometers. Applications of the theory include the use of ZR and RZ correlations within the SPAC methodology and the demonstration of increased robustness for ZR and RZ correlations in the presence of anisotropic noise distributions. A clear relationship exists between the time-domain correlations of ambient noise and the results of the SPAC method for the entire matrix of three-component correlation coefficients. Future topics to be explored include the extension of the 1-D and 2-D theories presented here to the 3-D case with P and S waves, inversion of the correlation spectra for the noise distribution, and the use of ZR and RZ correlations for imaging and monitoring the structure of the subsurface.

## ACKNOWLEDGMENTS

MH thanks Professor Paul Michaels of Boise State University for initially giving him a copy of the 1957 paper by Aki. This manuscript has benefitted from comments by Bernard Chouet, Lapo Boschi, and an anonymous reviewer. This work was partially funded under NSF award EAR-1142154 and USGS/ARRA award G10AC00016.

## REFERENCES

Aki, K., 1957. Space and time spectra of stationary stochastic waves, with special reference to microtremors, *Bull. Earthq. Res. Inst. Tokyo Univ.*, **25**, 415-457.

Aki, K., 1965. A note on the use of microseisms in determining the shallow structures of the Earth's crust, *Geophysics*, **30**, 665-666.

Aki, K. & Richards, P., 1980. *Quantitative Seismology*, 1st edn, W. H. Freeman and Company, San Francisco, CA.

Aki, K. & Richards, P., 2002. *Quantitative Seismology*, 2nd edn, W. H. Freeman and Company, San Francisco.

Campillo, M. & Paul, A., 2003. Long-range correlations in the diffuse seismic coda, *Science*, **299**, 547-549.

Chávez-García, F.J. & Luzón, F., 2005. On the correlation of seismic microtremors, *J. geophys. Res.*, **110**, B11313, doi:10.1029/2005JB003671.

Chávez-García, F.J., Rodríguez, M. & Stephenson, W.R., 2005. An alternative approach to the SPAC analysis of microtremors: exploiting stationarity of noise, *Bull. seism. Soc. Am.*, **95**, 277-293.

Cho, I. & Tada, T., 2006. A generic formulation for microtremor exploration methods using three-component records from a circular array, *Geophys. J. Int.*, **165**, 236-258.

Chouet, B., 1996. New methods and future trends in seismological volcano monitoring, in *Monitoring and Mitigation of Volcano Hazards*, pp. 23-97, eds Scarpa, R. & Tilling, R.I., Springer, New York, NY.

Chouet, B., De Luca, G., Milana, G., Dawson, P., Martini, M. & Scarpa, R., 1998. Shallow velocity structure of Stromboli Volcano, Italy, derived from small-aperture array measurements of strombolian tremor, *Bull. seism. Soc. Am.*, **88**, 653-666.

Chouet, B., Dawson, P. & Arciniega-Ceballos, A., 2005. Source mechanism of Vulcanian degassing at Popocatépetl Volcano, Mexico, determined from waveform inversion of very long period signals, *J. geophys. Res.*, **110**, B07301, doi:10.1029/2004JB003524.

Cox, H., 1973. Spatial correlation in arbitrary noise fields with application to ambient sea noise, *J. acoust. Soc. Am.*, **54**, 1289-1301.

Dixon, J.P., Stihler, S.D. & Power, J.A., 2008. *Catalog of Earthquake Hypocenters at Alaskan Volcanoes: January 1 through December 31, 2007*, U.S. Geol. Surv. Data Ser., No. 367, U.S. Geological Survey, Reston, VA.

Dorman, J. & Ewing, M., 1962. Numerical inversion of seismic surface wave dispersion data and crust-mantle structure in the New York-Pennsylvania area, *J. geophys. Res.*, **67**, 5227-5241.

Durand, S., Montagner, J.P., Roux, P., Brenguier, F., Nadeau, R.M. & Ricard, Y., 2011. Passive monitoring of anisotropy change associated with the Parkfield 2004 earthquake, *Geophys. Res. Lett.*, **38**, L13303, doi:10.1029/2011GL047875.

Ekström, G., Nettles, M. & Abers, G.A., 2003. Glacial earthquakes, *Science*, **302**(5645), 622-624.

Ekström, G., Abers, G.A. & Webb, S.C., 2009. Determination of surface-wave phase velocities across US Array from noise and Aki's spectral formulation, *Geophys. Res. Lett.*, **36**, doi:10.1029/2009GL039131.

Ferrazzini, V., Aki, K. & Chouet, B., 1991. Characteristics of seismic waves composing Hawaiian volcanic tremor and gas-piston events observed by a near-source array, *J. geophys. Res.*, **96**, 6199-6209.

Haney, M.M., 2009. Infrasonic ambient noise interferometry from correlations of microbaroms, *Geophys. Res. Lett.*, **36**, L19808, doi:10.1029/2009GL040179.

Haney, M.M., 2010. Location and mechanism of very long period tremor during the 2008 eruption of Okmok Volcano from interstation arrival times, *J. geophys. Res.*, **115**, B00B05, doi:10.1029/2010JB007440.

Haney, M.M., Nies, A., Masterlark, T., Needy, S. & Pedersen, R., 2011. Interpretation of Rayleigh-wave ellipticity observed with multicomponent passive seismic interferometry at Hekla Volcano, Iceland, *Leading Edge*, **30**, 526-531.

Harmon, N., Rychert, C. & Gerstoft, P., 2010. Distribution of noise sources for seismic interferometry, *Geophys. J. Int.*, **183**, 1470-1484.

Lin, F.-C., Moschetti, M.P. & Ritzwoller, M.H., 2008. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps, *Geophys. J. Int.*, **173**, 281-298.
Lobkis, O.I. & Weaver, R.L., 2001. On the emergence of the Greens function in the correlations of a diffuse field, *J. acoust. Soc. Am.*, **110**, 3011-3017.

Lu, Z., Wicks, C., Jr, Power, J.A. & Dzurisin, D., 2000. Ground deformation associated with the March 1996 earthquake swarm at Akutan volcano, Alaska, revealed by satellite radar interferometry, *J. geophys. Res.*, **105**(B9), 21 483–21 495.

Masterlark, T., Haney, M., Dickinson, H., Fournier, T. & Searcy, C., 2010. Rheologic and structural controls on the deformation of Okmok Volcano, Alaska: FEMs, InSAR, and ambient noise tomography, *J. geophys. Res.*, **115**, B02409, doi:10.1029/2009JB006324.

Muyzert, E., 2007a. Seabed property estimation from ambient-noise recordings. Part 1: compliance and Scholte wave phase-velocity measurements, *Geophysics*, **72**, U21–U26.

Muyzert, E., 2007b. Seabed property estimation from ambient-noise recordings. Part 2: Scholte-wave spectral-ratio inversion, *Geophysics*, **72**, U47–U53.

Nakahara, H., 2006. A systematic study of theoretical relations between spatial correlation and Green's function in one-, two-, and three-dimensional random scalar wavefields, *Geophys. J. Int.*, **167**, 1097–1105.

Okada, H., 2003. *The Microtremor Survey Method*, Geophys. Monogr. No. 12, Society of Exploration Geophysicists, Tulsa, OK.

Prieto, G.A., Lawrence, J.F. & Beroza, G.C., 2009. Anelastic Earth structure from the coherency of the ambient seismic field, *J. geophys. Res.*, **114**, B07303, doi:10.1029/2008JB006067.

Rhie, J. & Romanowicz, B., 2004. Excitation of Earth's continuous free oscillations by atmosphere-ocean-seafloor coupling, *Nature*, **431**(7008), 552–556.

Roux, P., 2009. Passive seismic imaging with directive ambient noise: application to surface waves and the San Andreas Fault in Parkfield, CA, *Geophys. J. Int.*, **179**, 367–373.

Roux, P., Wathelet, M. & Roueff, A., 2011. The San Andreas Fault revisited through seismic-noise and surface-wave tomography, *Geophys. Res. Lett.*, **38**, L13319, doi:10.1029/2011GL047811.

Saccorotti, G., Chouet, B. & Dawson, P., 2003. Shallow-velocity models at the Kilauea Volcano, Hawaii, determined from array analyses of tremor wavefields, *Geophys. J. Int.*, **152**, 633–648.

Snieder, R., 2001. *A Guided Tour of Mathematical Methods for the Physical Sciences*, Cambridge University Press, Cambridge.

Snieder, R., 2002. Scattering of surface waves, in *Scattering and Inverse Scattering in Pure and Applied Science*, pp. 562–577, eds Pike, R. & Sabatier, P., Academic Press, San Diego, CA.

Snieder, R., Fan, Y., Slob, E. & Wapenaar, K., 2010. Equipartitioning is not sufficient for Green's function extraction, *Earthq. Sci.*, **23**, 403–415.

Stehly, L., Campillo, M. & Shapiro, N., 2007. Travel time measurements from noise correlation: stability and detection of instrumental errors, *Geophys. J. Int.*, **171**, 223–230.

Tada, T. & Cho, I., 2007. Beyond the SPAC method: Exploiting the wealth of circular-array methods for microtremor exploration, *Bull. seism. Soc. Am.*, **97**, 2080–2095.

Tanimoto, T. & Alvizuri, C., 2006. Inversion of the HZ ratio of microseisms for S-wave velocity, *Geophys. J. Int.*, **165**, 323–335.

Teal, P.D., Abhayapala, T.D. & Kennedy, R.A., 2002. Spatial correlation for general distributions of scatterers, *IEEE Signal Process. Lett.*, **9**, 305–308.

Tsai, V.C., 2009. On establishing the accuracy of noise tomography travel-time measurements in a realistic medium, *Geophys. J. Int.*, **178**, 1555–1564.

Tsai, V.C. & Moschetti, M.P., 2010. An explicit relationship between time-domain noise correlation and spatial autocorrelation (SPAC) results, *Geophys. J. Int.*, **182**, 454–460.

Van Wijk, K., Mikesell, T.D., Schulte-Pelkum, V. & Stachnik, J., 2011. Estimating the Rayleigh-wave impulse response between seismic stations with the cross terms of the Green tensor, *Geophys. Res. Lett.*, **38**, L16301, doi:10.1029/2011GL047442.

Zobin, V.M., Plascencia, I., Reyes, G. & Navarro, C., 2009. The characteristics of seismic signals produced by lahars and pyroclastic flows: Volcán de Colima, México, *J. Volc. Geotherm. Res.*, **179**(1-2), 157–167.

## APPENDIX A: GENERAL FORM FOR ANISOTROPIC INCIDENCE

For Rayleigh waves, the expression for the correlation function for generally anisotropic incidence is the following:

$$ \phi_{ij} = P^R(\omega_0) M_{ij} \sum_{m=0}^{\infty} \epsilon_m i^m J_m \left(\frac{\omega_0}{c} r \right) \text{Re} [\gamma_m^{ij}], \quad (A1) $$

where $P^R(\omega_0)$ is the power spectrum of the Rayleigh waves, the indices i, j take on values of either Z, R or T, $M_{ij}$ is a multiplier given by

$$ M_{ZZ} = 1 $$

$$ M_{ZR} = M_{ZT} = Ri \operatorname{sgn}(\omega_0) $$

$$ M_{RZ} = M_{TZ} = -Ri \operatorname{sgn}(\omega_0) $$

$$ M_{RR} = M_{RT} = M_{TR} = M_{TT} = R^2 $$

(A2)

and the angular integrals are given by

$$ \gamma_m^{ZZ} = \frac{1}{2\pi} \int_0^{2\pi} p(\theta) \exp[-im(\theta - \psi)] d\theta $$

$$ \gamma_m^{ZR} = \gamma_m^{RZ} = \frac{1}{2\pi} \int_0^{2\pi} p(\theta) \cos(\theta - \psi) \exp[-im(\theta - \psi)] d\theta $$

$$ \gamma_m^{ZT} = \gamma_m^{TZ} = \frac{1}{2\pi} \int_0^{2\pi} p(\theta) \sin(\theta - \psi) \exp[-im(\theta - \psi)] d\theta $$

$$ \gamma_m^{RR} = \frac{1}{2\pi} \int_0^{2\pi} p(\theta) \cos^2(\theta - \psi) \exp[-im(\theta - \psi)] d\theta $$

$$ \gamma_m^{TT} = \frac{1}{2\pi} \int_0^{2\pi} p(\theta) \sin^2(\theta - \psi) \exp[-im(\theta - \psi)] d\theta $$

$$ \gamma_m^{RT} = \gamma_m^{TR} = \frac{1}{2\pi} \int_0^{2\pi} p(\theta) \sin(\theta - \psi) \cos(\theta - \psi) \exp[-im(\theta - \psi)] d\theta. \quad (A3) $$
In the above expressions, R is the ratio of the horizontal-to-vertical motion. Also, note that the factor isgn ($\omega_0$) does not appear for $M_{RR}$, $M_{RT}$, $M_{TR}$ and $M_{TT}$ since the radial and transverse motions are not elliptically polarized, as is the case for the vertical and either of the horizontal motions.

For Love waves, the vertical displacement is zero and the horizontal displacements are

$$u_r(r, \psi, t) = \sum_n \sum_m A_{nm} \sin(\theta_m - \psi) \exp[i\rho_n r \cos(\theta_m - \psi)] \sin(c\rho_n t) + \sum_n \sum_m \frac{B_{nm}}{c\rho_n} \sin(\theta_m - \psi) \exp[i\rho_n r \cos(\theta_m - \psi)] \cos(c\rho_n t), \quad (A4)$$

and

$$u_{\psi}(r, \psi, t) = \sum_{n} \sum_{m} A_{nm} \cos(\theta_{m} - \psi) \exp[i\rho_{n}r \cos(\theta_{m} - \psi)] \sin(c\rho_{n}t) + \sum_{n} \sum_{m} \frac{B_{nm}}{c\rho_{n}} \cos(\theta_{m} - \psi) \exp[i\rho_{n}r \cos(\theta_{m} - \psi)] \cos(c\rho_{n}t). \quad (A5)$$

Applying the analysis to the RR, RT, TR and TT correlations gives the following general form

$$\phi_{ij} = P^{L}(\omega_0) \sum_{m=0}^{\infty} \epsilon_m i^m J_m \left(\frac{\omega_0}{c} r \right) \text{Re} [\gamma_m^{ij}] \quad (A6)$$

where $P^L(\omega_0)$ is the power spectrum of the Love waves, the indices i, j take on values of either R or T, and the angular integrals are given by

$$\begin{align*}
\gamma_m^{RR} &= \frac{1}{2\pi} \int_0^{2\pi} p(\theta) \sin^2(\theta - \psi) \exp[-im(\theta - \psi)] d\theta \\
\gamma_m^{TT} &= \frac{1}{2\pi} \int_0^{2\pi} p(\theta) \cos^2(\theta - \psi) \exp[-im(\theta - \psi)] d\theta \\
\gamma_m^{RT} &= \gamma_m^{TR} = \frac{1}{2\pi} \int_0^{2\pi} p(\theta) \sin(\theta - \psi) \cos(\theta - \psi) \exp[-im(\theta - \psi)] d\theta. \tag{A7}
\end{align*}$$

## APPENDIX B: ANGLE INTEGRALS OVER AN INTERVAL

For incident energy that is constant ($p(\theta) = 1$) over an angular range of $2\Delta$ centred on the angle $\phi_0$, the angular integrals are given by the following expressions. The azimuth between the seismometers is denoted by $\psi$.

$$\gamma_m^{ZZ} = \frac{1}{2\Delta} \int_{\phi_0 - \Delta}^{\phi_0 + \Delta} \exp[-im(\theta - \psi)] d\theta = \exp[-im(\phi_0 - \psi)] \text{sinc}(m\Delta). \quad (B1)$$

$$\gamma_m^{ZR} = \frac{1}{2\Delta} \int_{\phi_0 - \Delta}^{\phi_0 + \Delta} \cos(\theta - \psi) \exp[-im(\theta - \psi)] d\theta = \frac{\exp[-i(m-1)(\phi_0 - \psi)] \text{sinc}[(m-1)\Delta]}{2} + \frac{\exp[-i(m+1)(\phi_0 - \psi)] \text{sinc}[(m+1)\Delta]}{2}. \quad (B2)$$

$$\gamma_m^{ZT} = \frac{1}{2\Delta} \int_{\phi_0 - \Delta}^{\phi_0 + \Delta} \sin(\theta - \psi) \exp[-im(\theta - \psi)] d\theta = -i \frac{\exp[-i(m-1)(\phi_0 - \psi)] \text{sinc}[(m-1)\Delta]}{2} + i \frac{\exp[-i(m+1)(\phi_0 - \psi)] \text{sinc}[(m+1)\Delta]}{2}. \quad (B3)$$

$$\gamma_m^{RT} = \frac{1}{2\Delta} \int_{\phi_0 - \Delta}^{\phi_0 + \Delta} \sin(\theta - \psi) \cos(\theta - \psi) \exp[-im(\theta - \psi)] d\theta = -i \frac{\exp[-i(m-2)(\phi_0 - \psi)] \text{sinc}[(m-2)\Delta]}{4} + i \frac{\exp[-i(m+2)(\phi_0 - \psi)] \text{sinc}[(m+2)\Delta]}{4}. \quad (B4)$$

$$\gamma_m^{RR} = \frac{1}{2\Delta} \int_{\phi_0 - \Delta}^{\phi_0 + \Delta} \cos^2(\theta - \psi) \exp[-im(\theta - \psi)] d\theta = \frac{\exp[-i(m-2)(\phi_0 - \psi)] \text{sinc}[(m-2)\Delta]}{4} + \frac{\exp[-im(\phi_0 - \psi)] \text{sinc}(m\Delta)}{2} + \frac{\exp[-i(m+2)(\phi_0 - \psi)] \text{sinc}[(m+2)\Delta]}{4}. \quad (B5)$$

$$\gamma_m^{TT} = \frac{1}{2\Delta} \int_{\phi_0 - \Delta}^{\phi_0 + \Delta} \sin^2(\theta - \psi) \exp[-im(\theta - \psi)] d\theta = -\frac{\exp[-i(m-2)(\phi_0 - \psi)] \text{sinc}[(m-2)\Delta]}{4} + \frac{\exp[-im(\phi_0 - \psi)] \text{sinc}(m\Delta)}{2} - \frac{\exp[-i(m+2)(\phi_0 - \psi)] \text{sinc}[(m+2)\Delta]}{4}. \quad (B6)$$
APPENDIX C: CONNECTION TO GREEN'S FUNCTION FOR $\omega < 0$ IN 2-D

For $\omega < 0$, we find that

$$ -R J_1\left(\frac{|\omega_0|}{c}r\right) = \frac{R}{2}\sqrt{\frac{2c}{\pi|\omega|r}}\left[\exp\left(-\frac{i|\omega|}{c}r\right) + i\exp\left(\frac{i|\omega|}{c}r\right)\right]\exp\left(-\frac{i\pi}{4}\right). \quad (C1) $$

For the Green's function, we take the following expression from the matrix in eq. (70)

$$ G_{ZR} = \frac{1}{4} \times \sqrt{\frac{2c}{\pi|\omega|r}} R \exp\left(-\frac{i|\omega|}{c}r - i\frac{3\pi}{4}\right) \quad (C2) $$

to find that

$$ -2(G_{ZR} - G_{ZR}^*) = -\frac{R}{2}\sqrt{\frac{2c}{\pi|\omega|r}}\left[\exp\left(\frac{i|\omega|}{c}r\right) - i\exp\left(-\frac{i|\omega|}{c}r\right)\right]\exp\left(-\frac{i\pi}{4}\right). \quad (C3) $$

Eqs (C1) and (C3) can be made equal to each other if eq. (C1) is multiplied by a factor of $i$. This factor is the Hilbert transform operator –isgn($\omega$) for the case we are currently considering, $\omega < 0$.