# Incoherent noise-induced distortions of Rayleigh wave ellipticity measurements obtained with three-component beamforming

Marc Wathelet

► To cite this version:

Marc Wathelet. Incoherent noise-induced distortions of Rayleigh wave ellipticity measurements obtained with three-component beamforming. Geophysical Journal International, 2024, 236 (3), pp.1804-1827. 10.1093/gji/ggae017. hal-04908185

HAL Id: hal-04908185

https://hal.science/hal-04908185v1

Submitted on 21 Mar 2025

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>

Distributed under a Creative Commons Attribution 4.0 International License
# Incoherent noise-induced distortions of Rayleigh wave ellipticity measurements obtained with three-component beamforming

Marc Wathelet

Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTERRE, 38000 Grenoble, France. E-mail: Marc.Wathelet@univ-grenoble-alpes.fr

Accepted 2024 January 10. Received 2023 December 5; in original form 2023 August 31

## SUMMARY

For site characterization, the elliptic particle motion of Rayleigh waves and its frequency dependence is a well-known property that aroused less interest than the frequency dependence of the phase velocity. More than 50 yr ago, ellipticity was already recognized as providing information independent from phase velocity, despite the difficulties inherent to its accurate and precise measurement. Several techniques were developed during the last two decades to extract the ellipticity curve from ambient vibration recordings, with a single three-component (3C) station, with pairs of 3C stations and more recently with 3C arrays. The latter has the advantage over the other approaches that the sign of the ellipticity can be retrieved. Moreover, higher order mode separation is possible under certain conditions. Nevertheless, Rayleigh Three-component BeamForming (RTBF) proposed by Wathelet et al. encounters difficulties in the presence of significant levels of incoherent noise when the true ellipticity is vanishing or when it has a high absolute value. In this work, the analytical expressions of the beam power for a single source wavefield are revised under more realistic assumptions for the incoherent noise azimuthal distribution. The proposed model also includes an asymmetric distribution of the incoherent noise between vertical and horizontal components, which was not the case in the original publication. Switching from ellipticity to angular ellipticity drastically simplifies the formalism. Moreover, it naturally leads to a new steering matrix (all-component ellipticity steering) which solves the limitation around zero and infinity observed for RTBF. Interestingly, the accuracy of the ellipticity is no longer influenced by the absolute level of incoherent noise but by the difference between the incoherent noise ratio on vertical and horizontal components. A method based on the second derivative of the beam power versus the radial wavenumber is finally proposed to experimentally measure the noise ratio difference, which allows experimental values to be corrected. The methodology is compared with classical vertical beamforming and RTBF for a synthetic case and three experimental data sets.

**Key words:** Fourier analysis; Time-series analysis; Surface waves and free oscillations; Wave propagation.

## 1 INTRODUCTION

Elliptic orbits are one of the Rayleigh surface wave properties discovered by Strutt (1885) for a half-space. Thereafter Love (1911) tried to analyse the effects of a superficial layer on Rayleigh waves for incompressible media. Because of the complexity of the equations, the interpretation was limited to a numerical example with a wavelength of the same order as the layer thickness. The amplitude ratio of vertical and horizontal displacements was found as nearly equivalent to the value found for a half-space. The frequency dependence of the ratio was not yet considered. Solving the wave equations to get the frequency dependence of the phase velocity for heterogeneous media (one or two superficial layers)

concentrated most of the efforts during the first half of the 20th century (Love 1911; Sezawa 1927; Stoneley 1934; Sezawa & Kanai 1938; Fu 1946; Scholte 1947). Thomson (1950) and Haskell (1953) first provided a matrix formulation that enabled a dispersion curve to be calculated in a reasonable time for an unlimited stack of layers. Many improvements were proposed during the following decades, for example to solve high frequency instabilities (Dunkin 1965).

The first exhaustive parametric study of the characteristics of Rayleigh wave dispersion for a uniform layer overlying a half-space can be found in Mooney & Bolt (1966). Namely, the ratio of horizontal to vertical surface displacements, readily available with phase velocity computation, is explicitly plotted as a function of period
and adimensional model parameters. Boore & Toksöz (1969) pro-
posed to use the phase velocity and the ellipticity for determining
earth structures, the latter contributing to independent information.
It was confirmed later by Hobiger et al. (2013) and Boaga et al.
(2013). Sexton et al. (1977) noted that only a few observations of
ellipticity are reported in the literature at that time. After observ-
ing repeatable and azimuthally independent ellipticity curves for
several earthquake events in the U.S. Midwest, they successfully
compared them with existing theoretical models, also showing the
influence of superficial layers.

The aforementioned references are mainly focused on earthquake
data. Microtremors became more and more popular at the end of the
20th century. In particular, the horizontal to vertical spectral ratios
(HVSR curve versus frequency), a simple but powerful technique,
focused a lot of attention after the work of Nakamura (1989). There
were intense—and there are probably still some—debates about the
origin of the HVSR peaks and the amplitude of the curve (Mol-
nar et al. 2018): theoretical 1-D transfer function of body waves
(Nakamura 2000), Rayleigh waves ellipticity (Lachet & Bard 1994;
Lermo & Chávez-García 1994; Fäh et al. 2001) or more recently
diffuse wavefields (Sánchez-Sesma et al. 2011). Fäh et al. (2001)
opened an alternative way by computing the ratio only for the por-
tions of signals where vertical–radial polarization dominates, that
is, where HVSR should tend to the ellipticity. They were followed
by Hobiger et al. (2009) who used a random decrement technique to
extract the ellipticity from single station recordings, taking advan-
tage of the constant theoretical phase difference between the radial
and the vertical component of a Rayleigh wave.

With a single three-component (3C) station, since there is no
information about the propagation azimuth, the decomposition be-
tween radial and transverse components is impossible without sup-
plementary assumptions. Hobiger et al. (2009) solves this issue
by maximizing the correlation between the vertical and the radial
component, which may encounter difficulties if either of the two
components vanishes and in the case of multiple Rayleigh modes.
Higher order modes are reported to explain the seasonality of el-
lipticity (Tanimoto et al. 2006). The same relationship between
polarization and azimuth is exploited by Takagi et al. (2018) to
obtain the dominant azimuth from single station processing. In the
case of a directive seismic-noise distribution, Roux (2009) produced
an ellipticity map for the Parkfield area (United States) thanks to
correlation tensors calculated from each pair of 3C stations. A pair
of stations contains supplemental information compared to a sin-
gle station: a more robust separation of Rayleigh and Love waves
might be expected. Nonetheless, 30 stations and days of signals were
necessary to achieve these results for a wavefield with a dominant
propagation direction.

Microtremors, with impinging waves produced by a diversity of
simultaneous sources, not necessarily with a dominant source nor
a uniform azimuthal distribution, have been successfully processed
by beamforming to obtain phase velocities versus frequency (La-
coss et al. 1969; Asten & Henstridge 1984; Garofalo et al. 2016a).
Most of these works were including the vertical component only.
The horizontal components, which record a mixture of Love and
Rayleigh waves, were generally considered as too complex com-
pared to the vertical components dominated by Rayleigh waves. Fäh
et al. (2008) proposed an efficient separation of Rayleigh and Love
waves with a high-resolution beamforming (Capon 1969) based on
a cross-spectral matrix constructed with the projection of the hori-
zontal components on the presumed direction and its perpendicular
for Rayleigh and Love, respectively. Poggi & Fäh (2010) combined
the vertical and the radial beam powers to produce experimental

unsigned ellipticity curves. The two types of polarization are also
separated with standard beamforming (Lacoss et al. 1969) and a
cross-spectral matrix made of the 3N raw components if N is the
number of sensors (Riahi et al. 2013). The sign and the amplitude of
the ellipticity were modelled through a set of 91 predefined polar-
ization states. Each state corresponds to a phase shift and a relative
amplitude modulation of the original signal components. This study
processed arrays made of about 80 3C sensors. Gal et al. (2016)
and Löer et al. (2018) also used 3C cross-spectral matrices to study
ocean induced microseisms and anisotropy, respectively. Elliptic-
ity curves are not shown nor analysed in these three publications.
MUSIQUE (Hobiger et al. 2012) and WaveDec (Maranò et al. 2012,
2017) are also two array-based methods to extract signed ellipticity
curves. Hobiger et al. (2021) provides a review of these techniques.

Using a low number of sensors (about 15) and high-resolution
beamforming, Wathelet et al. (2018) built a 2N cross-spectral matrix
with vertical and horizontal components, the latter being projected
on the presumed radial direction like in Fäh et al. (2008). By con-
sidering analytically a wavefield made of a single plane wave, they
suggested two ways to define the steering matrix involved in a 3C
beamforming, namely, $E'_{Rh}$ and $E'_{Rz}$ (following their notations). The
ellipticity estimators found for both steering matrices were found to
be biased by the level of incoherent noise. A third way of comput-
ing the ellipticity was proposed to remove the effect of incoherent
noise.

In the present work, the assumptions made about the incoherent
noise distribution in Wathelet et al. (2018) will be reconsidered. A
more realistic incoherent noise model including distinct noise levels
for vertical and horizontal components will be proposed. The analyt-
ical formulations of the beam power will be developed for a single
plane wave using 2N [Rayleigh Three-component BeamForming
(RTBF)] and 3N cross-spectral matrices [Rayleigh Direct Steering
(RDS)]. For both approaches, a formalism based on an angular el-
lipticity ($\xi$ = arctan e where e is the ellipticity in the usual sense)
will be introduced which will naturally gives rise to the concept
of All-component ellipticity steering (A- prefix to RTBF and RDS
acronyms). The classical nested grid search usually implemented for
2-D FK ($k_x$ and $k_y$) is not very efficient when an additional parame-
ter is deemed. Therefore, a fast algorithm based on the beam power
derivative versus ellipticity will be presented for searching the com-
bined wavenumber-ellipticity space. The analysis of the properties
of the proposed ellipticity estimators will show that the incoherent
noise level and its partition between horizontal and vertical compo-
nents must be accurately measured to retrieve bias-free ellipticity
curves. However, we will see that the maximization of the beam
power offers no way to determine these parameters. A technique
based on the second derivative of the beam power will be proposed
to assess the incoherent noise parameters and consequently to cor-
rect the ellipticity values. The various array-based approaches for
the retrieval of Rayleigh wave ellipticity curves mentioned herein-
above will be discussed and compared for one synthetic case and
three experimental sites.

## 2 DEFINITIONS

We assume that an array of *N* 3C sensors, located on a flat ground surface at $\vec{r}_i$, is recording a Rayleigh plane wave propagating at wavenumber $\vec{k}_r$. The wavenumber is in the horizontal plane with the direction $\theta_r$, counted positively from east to north axes. For the three axes *E*, *N* and *Z* (positive upwards), the Fourier spectra of the
signals observed at each sensor for the angular frequency $\omega$ are

$$ 
\begin{aligned}
X_E(\omega, \vec{k}_r) &= S_{R_h}(\omega) \cos \theta_r \mathbf{q}(\vec{k}_r) + \eta_E(\omega) \\
X_N(\omega, \vec{k}_r) &= S_{R_h}(\omega) \sin \theta_r \mathbf{q}(\vec{k}_r) + \eta_N(\omega) \\
X_Z(\omega, \vec{k}_r) &= S_{R_z}(\omega) \mathbf{q}(\vec{k}_r) + \eta_Z(\omega),
\end{aligned}
\quad (1) $$

where $S_{R_z}$ and $S_{R_h}$ are the vertical and the radial spectra of the incident plane wave, respectively; $\eta_E$, $\eta_N$ and $\eta_Z$ are the incoherent noise spectra recorded at each sensor for the three axes ($N \times 1$ matrices); $\mathbf{q}(\vec{k})$ is the phase shift matrix ($N \times 1$) corresponding to the propagation at wavenumber $\vec{k}$: $q_i(\vec{k}) = \exp(-j\vec{k} \cdot \vec{r}_i)$, where $j$ is the imaginary unit. For the sake of conciseness, $\mathbf{q}(\vec{k})$ and $\mathbf{q}(\vec{k}_r)$ are written $\mathbf{q}$ and $\mathbf{q}_r$, respectively. Bold face characters are matrices.

A beamforming technique requires the computation of the averaged cross-spectrum density matrix ($F$) over a set of $M$ time blocks (Capon 1969). If only the vertical component of the sensors is considered, we have

$$ 
\begin{aligned}
F_z(\omega, \vec{k}_r) &= X_Z(\omega, \vec{k}_r) X_Z^*(\omega, \vec{k}_r) \\
&= |S_{R_z}(\omega)|^2 \mathbf{q}(\vec{k}_r) \mathbf{q}^*(\vec{k}_r) + |S_{N_z}(\omega)|^2 \mathbf{I},
\end{aligned}
\quad (2) $$

where $S_{N_z}$ is the average incoherent noise spectrum on the vertical component, assumed to be the same for all sensors; $\mathbf{I}$ is the identity matrix. In Capon (1969), a ratio $R$ is defined to measure the relative amount of incoherent noise compared to the coherent wave ($R = |S_{N_z}|^2 / |S_{R_z}|^2$). It allows factoring out $|S_{R_z}|^2$ in eq. (2) which leads to a simple expression that depends only upon $\mathbf{q}_r$ and $R$. In the context of a 3C wavefield, it is shown in the next section that a slightly different definition of $R$ must be adopted.

## 2.1 A 3C incoherent noise model

In Wathelet et al. (2018), the ratio of the incoherent noise over the coherent signal $R$ for a noisy wavefield dominated by a Rayleigh plane wave is defined as

$$ R = \frac{|S_{N_h}|^2}{\cos^2(\theta - \theta_r) |S_{R_h}|^2}, \quad (3) $$

where $S_{N_h}$ is the noise spectrum observed for any arbitrary horizontal direction $\theta$. It is assumed to be the same in all directions and for all sensors (Wathelet et al. 2018). If $|S_{N_h}|^2$ is constant for all $\theta$, $R$ is thus a function of $\theta$ because of the cosine. When $\theta$ is equal to $\theta_r$, $R$ has a minimum equal to $|S_{N_h}|^2 / |S_{R_h}|^2$. On both sides of this minimum, $R$ tends to infinity when the angular difference tends to $\pm 90^\circ$. In Wathelet et al. (2018), the global maximum of the beam power would not be at $\vec{k} = \vec{k}_r$ if $R$ had such a dependence on $\theta$. Hence, $R$ is necessarily considered as independent of $\theta$ in Wathelet et al. (2018) which implicitly requires a dependence of $S_{N_h}$ on $\theta$. They assume that there is no incoherent noise on the directions perpendicular to the direction of propagation $\theta_r$. This is a very particular case, certainly not usually encountered for experimental wavefields. Therefore, formulating the assumptions of the problem in a different way is necessary. $R$ cannot be expressed for a single component like in eq. (3) because the coherent signal in any arbitrary direction $\theta$ is not constant, or may even vanish for particular combinations between the direction of propagation and the polarization. Taking into account the contributions of all components, the total energy of a Rayleigh harmonic plane wave remains constant (without scattering or intrinsic absorption) even if the Rayleigh signal may be null

for some particular directions or components. We can write that the incoherent noise ratio $R$ is

$$ R = \frac{|S_{N_e}|^2 + |S_{N_h}|^2 + |S_{N_z}|^2}{|S_{R_z}|^2 + |S_{R_h}|^2} = \frac{|S_{N_e}|^2 + |S_{N_h}|^2 + |S_{N_z}|^2}{|S_R|^2}, \quad (4) $$

where $|S_R|^2$ is the total spectrum power of the coherent wave; $S_{N_e}$ and $S_{N_h}$ are the incoherent noise spectra for the east and the north components, respectively. We can also define the same ratio restricted to the vertical or to the horizontal components:

$$ R_e = \frac{|S_{N_e}|^2}{|S_R|^2}, \quad R_n = \frac{|S_{N_h}|^2}{|S_R|^2}, \quad R_z = \frac{|S_{N_z}|^2}{|S_R|^2}. \quad (5) $$

$R$ is the sum of these three ratios. If the incoherent noise wavefield is equally distributed in all horizontal directions, the incoherent noise ratios $R_e$ and $R_n$ have the same value $R_h$, which is assumed in following sections.

## 2.2 Angular ellipticity

A Rayleigh wave is characterized by a particle motion in the vertical and in the radial directions. In an elastic medium, the phase difference between the two components is exactly $\pm 90^\circ$, which induces a rotation either retrograde or prograde according to the sign of the difference. For non-elastic cases, the absolute value of the phase difference is slightly changed which produced a tilt of the axes of the elliptical particle motion (Boore & Toksöz 1969). Therefore, the relationship between $S_{R_z}$ and $S_{R_h}$ is controlled by a complex ellipticity defined by (Wathelet et al. 2018, eq. 9):

$$ \frac{S_{R_h}}{S_{R_z}} = -e_r(\tau_r + j), \quad (6) $$

where $e_r$ is the spectral amplitude ratio (or the ellipticity) of the Rayleigh wave and $\tau_r$ is the ratio of the real over the imaginary parts of the ellipticity. Maranò et al. (2012) suggested the ellipticity angle to formulate the particle displacement generated by a Rayleigh wave. In a similar way, we can use an angular ellipticity to link the vertical spectrum $S_{R_z}(\omega)$ and the horizontal spectrum $S_{R_h}(\omega)$ to the global Rayleigh spectrum $S_R(\omega)$:

$$ 
\begin{aligned}
S_{R_z}(\omega) &= t_r(\omega)S_R(\omega) \cos \xi_r(\omega) \\
S_{R_h}(\omega) &= -t_r(\omega)(\tau_r(\omega) + j)S_R(\omega) \sin \xi_r(\omega),
\end{aligned}
\quad (7) $$

where

$$ 
\begin{aligned}
\xi_r(\omega) &= \arctan(e_r(\omega)) \\
t_r(\omega) &= \frac{1}{\sqrt{\cos^2 \xi_r(\omega) + (1 + \tau_r^2(\omega)) \sin^2 \xi_r(\omega)}}.
\end{aligned}
\quad (8) $$

$\xi_r(\omega)$ is the angular ellipticity; $t_r(\omega)$ is a factor that ensures that the sum of the vertical and the radial powers is equal to the total power. If $\tau_r^2$ is small compared to 1, $t_r$ becomes 1. The dependency to $\omega$ is dropped in the following sections for the sake of conciseness.

# 3 METHODS

## 3.1 Vertical beamforming (VBF)

Considering the angular ellipticity defined in eq. (8) and according to the definition of the incoherent noise ratios in eqs (4) and (5), the cross-spectrum density matrix for the vertical component in eq. (2) becomes

$$ F_z = |S_R|^2 (t_r^2 \cos^2 \xi_r \mathbf{q}_r \mathbf{q}_r^* + R_z \mathbf{I})^\dagger. \quad (9) $$
Symbol $^\dagger$ after an equation indicates that the steps used to obtain this equation are detailed in supplementary materials ('detailed equation steps'). Because $q_r^*q_r = N$, the Sherman-Morrison formula (Press et al. 2007) can be transformed to

$$(a q_r q_r^* + b I)^{-1} = \frac{1}{b} \left[ I - \frac{a}{Na + b} q_r q_r^* \right]^{\dagger}, \quad (10)$$

where a and b are real. Hence, an analytical form of the inverse of $F_z$ can be calculated.

$$F_z^{-1} = \frac{1}{|S_R|^2 R_z} \left[ I - \frac{t_r^2 \cos^2 \xi_r}{N t_r^2 \cos^2 \xi_r + R_z} q_r q_r^* \right] \quad (11)$$

The high-resolution beam power is defined like in Capon (1969):

$$P_z(\vec{k}) = \frac{1}{q^*(\vec{k}) F_z^{-1} q(\vec{k})}, \quad (12)$$

where $q(\vec{k})$ is the steering matrix for the presumed wavenumber $\vec{k}$. From the definition of $F_z$ and its inverse, we have

$$P_z(\vec{k}) = |S_R|^2 \frac{R_z}{N} \left[ 1 - \frac{|B(\vec{k} - \vec{k}_r)|^2}{1 + \frac{R_z}{N t_r^2 \cos^2 \xi_r}} \right]^{-1} \dagger. \quad (13)$$

$|B(\vec{k})|^2$ is the beamforming array response pattern normalized by the number of sensors. It depends only upon the array geometry $\vec{r}_i$. When there is a strong impedance contrast, the ellipticity is infinite at the resonance frequency (Malischewsky & Scherbaum 2004), $\xi_r$ tends to 90° and the beam power loses its dependency to the presumed wavenumber. When the Rayleigh horizontal component vanishes ($\xi_r = 0°$) and in the elastic case ($t_r = 1$), $P_z$ is the same as the beam power found by Capon (1969). If $\vec{k}$ is equal to $\vec{k}_r$, we have

$$P_z(\vec{k}_r) = |S_R|^2 \left( t_r^2 \cos^2 \xi_r + \frac{R_z}{N} \right). \quad (14)$$

## 3.2 Radial beamforming (RBF)

In a similar way, the cross-spectrum density matrix for the radial component is

$$F_h = |S_R|^2 (t_r^2 |\mu_\theta|^2 \sin^2 \xi_r q_r q_r^* + R_h I)^\dagger, \quad (15)$$

where

$$\mu_\theta = \cos(\theta - \theta_r)(\tau_r + j). \quad (16)$$

$F_h$ is the radial component computed from the horizontal components projected on the direction $\theta$ of the presumed wavenumber $\vec{k}$. From eqs (10) and (12) the beam power for the radial component is

$$P_h(\vec{k}) = |S_R|^2 \frac{R_h}{N} \left[ 1 - \frac{|B(\vec{k} - \vec{k}_r)|^2}{1 + \frac{R_h}{N t_r^2 |\mu_\theta|^2 \sin^2 \xi_r}} \right]^{-1} \dagger. \quad (17)$$

If $\vec{k}$ is equal to $\vec{k}_r$, it becomes

$$P_h(\vec{k}_r) = |S_R|^2 \left( t_r^2 (1 + t_r^2) \sin^2 \xi_r + \frac{R_h}{N} \right). \quad (18)$$

If we assume that $t_r^2$ is sufficiently small compared to 1, the ratio of $P_h$ over $P_z$ no longer provides the true ellipticity as suggested by Poggi & Fäh (2010).

$$e_p = \pm \sqrt{\frac{P_h(\vec{k}_r)}{P_z(\vec{k}_r)}} = \pm \sqrt{\frac{\sin^2 \xi_r + \frac{R_h}{N}}{\cos^2 \xi_r + \frac{R_z}{N}}} \quad (19)$$

For noise-free data sets, the ellipticity can be effectively retrieved by calculating the square root of the beam power ratio. However for experimental cases, the incoherent noise ratios are rarely negligible, particularly when either the sine or the cosine is small, that is, when $\xi_r$ tends to 0 or $\pm 90°$, respectively. For very noisy cases, this method can potentially return $\sqrt{\frac{R_h}{R_z}}$ instead of the true ellipticity. For example, if $R_z/N = 1$, $R_h/N = 1$ and $\xi_r$ has a usual value for a half-space (from $-35°$ to $-29°$, see Section 1 in Supporting Information), the observed angular ellipticity ($\xi_p = \arctan e_p$) varies from $\pm 39.9°$ to $\pm 41.7°$. If there is more incoherent noise on the horizontal components (e.g. due to sensor orientation errors), the error is larger ($\pm 48.4°$ to $\pm 49.7°$ if $R_h/N = 2$).

## 3.3 Rayleigh three-component beamforming (RTBF)

The RTBF method proposed by Wathelet et al. (2018) is revisited in this section with the more realistic definition of the incoherent noise ratio proposed in eqs (4) and (5). With the angular ellipticity formalism, the cross-spectral density matrix for RTBF is transformed to

$$F_2(\theta) = \begin{pmatrix} F_h & F_x \\ F_x^* & F_z \end{pmatrix}, \quad (20)$$

where

$$F_x = -|S_R|^2 t_r^2 \mu_\theta \sin \xi_r \cos \xi_r q_r q_r^* \dagger. \quad (21)$$

$F_2(\theta)$ is a two-component cross-spectrum matrix where horizontal components are projected on the radial direction $\theta$, that is, the direction of the presumed wavenumber $\vec{k}$; $F_z$ and $F_h$ are defined in eqs (9) and (15), respectively. Even if $\xi_r$ is null or $\pm 90°$, $F_2$ remains finite and invertible, which is not the case in Wathelet et al. (2018). According to Henderson & Searle (1981), a block-wise inverse of $F_2$ is

$$F_2^{-1} = \begin{pmatrix} F_h^{-1} + F_h^{-1} F_x A^{-1} F_x^* F_h^{-1} - F_h^{-1} F_x A^{-1} \\ -A^{-1} F_x^* F_h^{-1} A^{-1} \end{pmatrix}, \quad (22)$$

where

$$A = F_z - F_x^* F_h^{-1} F_x. \quad (23)$$

With the identity of eq. (10), we can calculate the inverse of $F_h$ and $A$, and the other blocks of $F_2^{-1}$,

$$F_2^{-1} = \frac{1}{R_h \lambda_\theta N |S_R|^2} \begin{pmatrix} \lambda_\theta N I - |\mu_\theta|^2 \sin^2 \xi_r q_r q_r^* & \sigma \mu_\theta \cos \xi_r \sin \xi_r q_r q_r^* \\ \sigma \mu_\theta^* \cos \xi_r \sin \xi_r q_r q_r^* & \sigma [\lambda_\theta N I - \sigma \cos^2 \xi_r q_r q_r^*] \end{pmatrix}^\dagger, \quad (24)$$

where

$$\lambda_\theta = \sigma \cos^2 \xi_r + |\mu_\theta|^2 \sin^2 \xi_r + \frac{R_h}{N t_r^2}, \quad (25) \\ \sigma = \frac{R_h}{R_z}.$$

The steering matrices $E_{2,h}$ and $E_{2,z}$ defined by Wathelet et al. (2018) can be reformulated in terms of a presumed angular ellipticity $\xi$:

$$E_{2,h} = \begin{pmatrix} -j \tan \xi q \\ q \end{pmatrix}, \quad E_{2,z} = \begin{pmatrix} -j q \\ \cot \xi q \end{pmatrix}. \quad (26)$$

The place of $j$ for $E_{2,z}$ is not the same as in Wathelet et al. (2018) but without any influence over a high resolution power (eq. 12). We can note that $E_{2,h} = \tan \xi E_{2,z}$. The angular frequency formulation
suggests a third option obtained by multiplying $E_{2,h}$ by the cosine of $\xi$, that we can call the all-component ellipticity steering matrix (RTBF becomes ARTBF):

$$E_{2,a} = \begin{pmatrix} -j \sin \xi \mathbf{q} \\ \cos \xi \mathbf{q} \end{pmatrix}. \qquad (27)$$

These definitions that may sound arbitrary are actually particular cases of an ideal steering matrix discussed in Section 2 of Supporting Information.

Defining the high resolution beam powers $P_{2,a}$, $P_{2,h}$ and $P_{2,z}$ in the same way as $P_z$ in eq. (12), we have

$$\begin{align}
P_{2,a}(\vec{k}, \xi) &= |S_R|^2 \frac{R_h}{N} \left[ \sin^2 \xi + \sigma \cos^2 \xi - \left( \left( \cos(\theta - \theta_r) \sin \xi_r \sin \xi + \sigma \cos \xi_r \cos \xi \right)^2 \right. \right. \nonumber \\
&\qquad \left. \left. + \cos^2 (\theta - \theta_r) \tau_r^2 \sin^2 \xi_r \sin^2 \xi \right) \frac{|B(\vec{k} - \vec{k}_r)|^2}{\lambda_{\theta}} \right]^{-1} \dagger, \nonumber \\
P_{2,h}(\vec{k}, \xi) &= \cos^2 \xi P_{2,a}, \nonumber \\
P_{2,z}(\vec{k}, \xi) &= \sin^2 \xi P_{2,a}. \tag{28}
\end{align}$$

If $\tau_r^2$ is sufficiently small and if $\vec{k}$ is equal to $\vec{k}_r$ ($\theta = \theta_r$), we have

$$P_{2,a}(\vec{k}_r, \xi) = |S_R|^2 \frac{R_h}{N} \left[ \sin^2 \xi + \sigma \cos^2 \xi - \frac{(\sin \xi_r \sin \xi + \sigma \cos \xi_r \cos \xi)^2}{\lambda_{\theta}} \right]^{-1} . \quad (29)$$

Additionally, if the incoherent noise is at the same level for the vertical and the horizontal components, the beam power reduces to

$$P_{2,a}(\vec{k}_r, \xi) = |S_R|^2 \frac{R_h}{N} \left[ 1 - \frac{\cos^2(\xi - \xi_r)}{1 + \frac{R_h}{N}} \right]^{-1}. \quad (30)$$

And if $\xi$ is equal to $\xi_r$, we have

$$P_{2,a}(\vec{k}_r, \xi_r) = |S_R|^2 \left(1 + \frac{R_h}{N}\right). \quad (31)$$

Wathelet et al. (2018) concluded that the 3C beamforming better attenuates the incoherent noise than a single component high-resolution FK. This is no longer noted when taking into account more realistic assumptions about the incoherent noise. R was always over 2N while being over N for the single component beamforming. In this work, $R_h$ and $R_z$ always appear over N for both single and 3C processing. However, we note that the peak amplitude of $P_{2,a}$ is independent of the ellipticity $\xi_r$ while the peak amplitude of $P_z$ and $P_h$ (eqs 14 and 18) is controlled by the ellipticity. The square trigonometric functions being always less or equal to 1, the influence of the incoherent noise is lower or equal for $P_{2,a}$ compared to $P_z$ and $P_h$.

From the analytical expression of the first two derivatives of the beam power versus $\xi$, we can assert that $P_{2,a}(\vec{k}_r, \xi)$, $P_{2,h}(\vec{k}_r, \xi)$ and $P_{2,z}(\vec{k}_r, \xi)$ have a single maximum at

$$\hat{\xi}_a = \frac{1}{2} \operatorname{atan2} \left( \sin 2\xi_r, \cos 2\xi_r + \delta \right)^{\dagger}, \quad (32)$$

$$\hat{\xi}_h = \arctan \left( \tan \xi_r \frac{1}{1 + \frac{R_z}{N \cos^2 \xi_r}} \right)^{\dagger} \text{ and} \quad (33)$$

$$\hat{\xi}_z = \arctan \left( \tan \xi_r \left( 1 + \frac{R_h}{N \sin^2 \xi_r} \right) \right)^{\dagger}, \quad (34)$$

respectively. $\delta$ is the difference between the incoherent noise on the vertical and on the horizontal components: $\frac{R_z}{N} - \frac{R_h}{N}$. atan2 has the

<figure><img src="image_2.png" /><figcaption>Figure 1. Distortions induced by the all-component ellipticity steering matrix ($E_{2,a}$) for a single source noisy wavefield with unbalanced incoherent noise ratios between horizontal and vertical components: the relationship between the estimated and the true ellipticity for $R_z > R_h$ (red lines) and $R_z < R_h$ (blue lines). The absolute difference $|\delta|$ equals to 0.1 (plain lines), 0.2 (dashed lines), 0.5 (dotted lines), 1 (dashed-dotted lines), 2 (dotted lines with crosses), 5 (dashed lines with crosses) and 10 (plain lines with crosses). The wavenumber is assumed to be perfectly resolved ($\vec{k} = \vec{k}_r$). The real part of the ellipticity is supposed to be null ($\tau_r = 0$).</figcaption></figure>

advantage over arctan that it returns an angle from -180° to 180°, hence keeping the same sign for $\xi_r$ and $\hat{\xi}_a$, which is a condition to get a negative concavity of the beam power. $\hat{\xi}_a$ is equal to the true ellipticity $\xi_r$ for any level of incoherent noise if $\delta$ is null, which is consistent with eq. (30). Fig. 1 summarizes the influence of the incoherent noise difference between the vertical and the horizontal components. If the incoherent noise is greater on the horizontal components than on the vertical components ($\delta < 0$), the absolute value of the ellipticity is always overestimated. Underestimation is always observed if $\delta > 0$. If the incoherent noise difference $|\delta|$ reaches a critical threshold of 1, the behaviour is unstable: $\hat{\xi}$ tends to ±90° for all ellipticities if $\delta < 0$ (blue curves) and it tends to 0° in the opposite case (red curves).

As already noted by Wathelet et al. (2018) the estimation is biased by the incoherent noise in opposite ways for $P_{2,h}$ and $P_{2,z}$. The more realistic definition of R in eq. (4) adds non-linearities to the ellipticity response compared to Wathelet et al. (2018), the deviations being themselves functions of $\xi_r$. Moreover, splitting of R into $R_z$ and $R_h$ shows that the two ellipticity estimations are affected by the incoherent noise of only one component, that is, $R_z$ or $R_h$ but not both. In Figs 2(a) and (b), the beam power is plotted as a function of the presumed ellipticity for $P_{2,h}$ and $P_{2,z}$, respectively. Various true angular ellipticities ($\xi_r$) are considered from -90° to 60°. The amplitude value of the peak is not constant for all true ellipticity values: for a constant level of the incoherent noise, $P_{2,h}$ is higher for an ellipticity with a small absolute value and it decreases when the ellipticity approaches ±90°. The contrary is observed for $P_{2,z}$ which vanishes for a null true ellipticity. The positions of the peaks are also clearly shifted compared to the true ellipticity marked by the coloured dots. For $P_{2,h}$, a variation of $R_h$ modifies the shape of the beam power function but its peak remains at the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** RTBF beam powers versus the presumed ellipticity for a single source noisy wavefield assuming that the wavenumber is perfectly resolved ($\bar{k} = \bar{k}_r$): (a) $P_{2,h}$, (b) $P_{2,z}$ and (c) $P_{2,s}$. The beam powers are calculated for various true angular ellipticities: −90 (blue lines), −60°, −45°, −30°, 0°, +30°, +45° and +60° (red lines). The positions of the true ellipticities are highlighted by the coloured dots at the top of the panels. The relative incoherent noise ratios $\frac{R_z}{N}$ and $\frac{R_h}{N}$ are equal to 0.5 for the continuous lines. The dashed lines are calculated for $\frac{R_h}{N} = 1$ and $\frac{R_z}{N} = 1$ in panels (a) and (b), respectively. The real part of the ellipticity is supposed to be null ($\tau_r = 0$).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Distortions of RTBF method for a single source noisy wavefield: (a) $P_{2,h}$, (b) $P_{2,z}$ and (c) $P_{2,s}$. The relationship between the estimated and the true ellipticity is shown for the relative incoherent noise ratios ($\frac{R_z}{N}$ and $\frac{R_h}{N}$) at values: 0.1, 1 and 10. In panel (c), $\frac{R_z}{N}$ and $\frac{R_h}{N}$ are identified by the line style and the colours (blue for $\frac{R_h}{N} < \frac{R_z}{N}$, red for $\frac{R_h}{N} > \frac{R_z}{N}$), respectively. For the grey curves, they are equal. The darker is the curve, the larger is the difference between $\frac{R_z}{N}$ and $\frac{R_h}{N}$. The wavenumber is assumed to be perfectly resolved ($\bar{k} = \bar{k}_r$). The real part of the ellipticity is supposed to be null ($\tau_r = 0$).</figcaption></figure>

same position. The same effect is observed on $P_{2,z}$ for a variation of $R_z$. The influence of the incoherent noise levels on the peak position are quantified in Figs 3(a) and (b): the non-linear effects are more pronounced for true ellipticities closer to ±90° for $P_{2,h}$ and closer to 0° for $P_{2,z}$. Therefore, the necessary condition to measure reliable ellipticities is a low level of incoherent noise on both components, which is not encountered for usual experimental environments. $P_{2,h}$ and $P_{2,z}$ must be associated because they are reliable on distinct true ellipticity ranges.

To remove the noise bias, Wathelet et al. (2018) proposed a combined beam power (renamed here $P_{2,s}$) made of the product of $P_{2,h}$ and $P_{2,z}$. Under the same conditions as eq. (29), the beam power has two maxima located symmetrically around zero at

$$ \hat{\xi}_s = \arctan \left( \pm \sqrt{\frac{\sin^2 \xi_r + \frac{R_h}{N}}{\cos^2 \xi_r + \frac{R_z}{N}}} \right)^{\dagger}. \quad (35) $$

The sign of $\xi$ has an influence on $P_{2,s}$ only because of the second term of the square brackets in eq. (29), namely, the product of sines. The product of cosines is always positive because angular ellipticities are by definition between −90° and 90°. When the

signs of $\xi$ and $\xi_r$ are the same, the term value is greater than when the signs differ. Therefore, the highest maximum is reached when $\hat{\xi}_s$ and $\xi_r$ have the same signs, confirming the observations in Wathelet et al. (2018) even if the incoherent noise context is modified. The level of incoherent noise but also its distribution between components control the amplitude difference between the two maxima. For instance, a large $\sigma$ lessen the sign influence and tends to an equivalence of the two maxima. Keeping $\sigma$ constant while increasing $R_h$ has a similar effect. The plots described for $P_{2,h}$ and $P_{2,z}$ are also available for $P_{2,s}$ in Figs 2(c) and 3(c). The bias is partially reduced but it is not removed, whereas Wathelet et al. (2018) claimed its removal with $P_{2,s}$. We can observe that for noisy situations the experimental ellipticities tend to fall around ±45° instead of returning the correct value. It can also be directly observed in eq. (35) where the ratio tends to 1 if the incoherent noise ratios are large compared to 1. Nevertheless, if a sufficiently low level of incoherent noise can be ensured, the results remain valid. If $\hat{\xi}_h$ and $\hat{\xi}_z$ are independent of $R_h$ and $R_z$, respectively, the advantage is lost for $\hat{\xi}_s$.

As far as ellipticity determination is concerned, the performances of $P_{2,a}$ and $P_{2,s}$ are compared for $\xi_r > 0$ and for $\delta < 0$ in Fig. 4.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Influence of the incoherent noise ratio difference on the vertical and the horizontal components ($\delta = R_z/N - R_h/N$) for $P_{2,a}$, $P_{2,s}$, $P_{2,h}$ and $P_{2,z}$: the difference between the estimated and the true ellipticities ($\hat{\xi}_? - \hat{\xi}_r$) is plotted versus the true ellipticity $\xi_r$ for various $\delta$ values: -0.1 (plain lines), -0.2 (dashed lines), -0.5 (dotted lines) and -1 (dashed-dotted lines). $\hat{\xi}_?$ is either $\hat{\xi}_a$ (black lines), $\hat{\xi}_s$ (blue lines without symbol for $R_z/N = 0.1$ and red lines without symbol for $R_z/N = 1$), $\hat{\xi}_h$ (plain lines with crosses) or $\hat{\xi}_z$ (plain lines with circles). $\hat{\xi}_?$ stands for $\hat{\xi}_a$, $\hat{\xi}_s$, $\hat{\xi}_h$ or $\hat{\xi}_z$. For $\hat{\xi}_h$, the blue line is calculated at $R_z/N = 0.1$ and the red line at $R_z/N = 1$. For $\hat{\xi}_z$, the blue line is calculated at $R_h/N = 0.1$ and the red line at $R_h/N = 1$. The line style does not apply for $\hat{\xi}_h$ and $\hat{\xi}_z$ because they are not influenced by $\delta$. The wavenumber is assumed to be perfectly resolved ($\bar{k} = \bar{k}_r$). The real part of the ellipticity is supposed to be null ($\tau_r = 0$).</figcaption></figure>

$P_{2,a}$ and $P_{2,s}$ with the same line style have the same $\delta$. For example, the three dashed lines correspond to the same incoherent noise difference between components ($\delta = -0.2$). $P_{2,a}$ is not influenced at all by the absolute noise level while it is the case for $P_{2,s}$: the comparison of the shape of the blue and the red curves gives an idea of the dependence of $P_{2,s}$ on the absolute noise level. If $|\delta|$ does not exceed 0.2, $\hat{\xi}_a$ deviates by a maximum of 6° around 45° and it is obviously the best choice. For $|\delta|$ between 0.2 and 0.5, and if $R_z/N$ is small enough, the smallest errors are observed for $\hat{\xi}_s$ in a broad range of true ellipticities around 45°. If $|\xi_r|$ is lower than 30°, $\hat{\xi}_a$ is preferred. For $|\delta|$ above 0.5, errors of more than 45° can be expected for all methods except around 0°, 45° and 90° for $\hat{\xi}_h$, $\hat{\xi}_s$ and $\hat{\xi}_z$, respectively. Under usually noisy environments all ellipticity estimators are biased. The best estimator depends upon the true ellipticity value which is precisely the target of the whole process.

Wathelet et al. (2018) proposed an experimental way to calculate the incoherent noise ratio assuming that it is the same on all components. As a consequence, noisy results could be filtered out which partially mitigates the imperfections of the method. Under a more realistic noise model, this calculation must be adjusted. From eqs (32) to (35), we have

$$ \frac{2}{\cot \hat{\xi}_h - \tan \hat{\xi}_z} = \frac{\sin 2\hat{\xi}_r}{\cos 2\hat{\xi}_r + \frac{R_z}{N} - \frac{R_h}{N}} = \tan 2\hat{\xi}_a^\dagger, \quad (36) $$

$$ \tan \hat{\xi}_h \tan \hat{\xi}_z = \tan^2 \hat{\xi}_s. \quad (37) $$

$\hat{\xi}_a$ and $\hat{\xi}_s$ are thus two ways to approach the true ellipticity that are directly related to $\hat{\xi}_h$ and $\hat{\xi}_z$. They are not independent quantities. From eqs (33) to (34), we can also deduce the following inequalities:

$$ |\hat{\xi}_h| \le |\xi_r| \le |\hat{\xi}_z|. \quad (38) $$

They are also valid if $\xi_r$ is replaced by $\hat{\xi}_a$ or $\hat{\xi}_s$ but the order of these three quantities depends upon $\xi_r$, $R_h/N$ and $R_z/N$ in a complex way. If $\delta$ is close to 0, $\hat{\xi}_a$ is a good approximation of $\xi_r$. From eqs (33) and (34), the incoherent noise level can be estimated with

$$ \frac{R_z}{N} = \frac{R_h}{N} = \frac{1}{2} \left[ \left( \cot \hat{\xi}_h + \tan \hat{\xi}_z \right) \frac{\sin 2\hat{\xi}_a}{2} - 1 \right]^\dagger. \quad (39) $$

The determination of the incoherent noise under this assumption does not bring substantial information about any potential bias of $\hat{\xi}_a$: theoretically, $\hat{\xi}_a$ should perfectly match the true ellipticity. Eq. (32) shows that the unique factor that control the accuracy of $\hat{\xi}_a$ is $\delta$. Hence a perceptive analysis requires an experimental assessment of $\delta$. There are thus three unknowns ($\xi_r$, $R_h$ and $R_z$) but only two independent equations. As a consequence, a reliable ellipticity under noisy conditions is impossible to attain by the maximization of a 3C beam power projected on the radial direction like $P_{2,a}$, $P_{2,h}$, $P_{2,z}$ and $P_{2,s}$. In the next section, we will explore the capabilities of a 3C cross-spectral matrix without projections like $F_2$: the radial projection is moved to the steering matrix.

### 3.4 All-component Rayleigh direct steering (ARDS)

A 3C cross-spectral density matrix $F_3$ is built from the raw components $X_E$, $X_N$ and $X_Z$, in the same way as Wagner (1996) or Riahi et al. (2013).

$$ F_3 = \frac{1}{M} \sum_{m=1}^{M} \left[ \begin{pmatrix} X_{E,m} \\ X_{N,m} \\ X_{Z,m} \end{pmatrix} (X_{E,m}^* X_{N,m}^* X_{Z,m}^*) \right] \quad (40) $$

From the definitions in eq. (1) for block $m$, the relations between $S_{R_z}$, $S_{R_h}$ and $S_R$ in eq. (8) and assuming that $M$ is sufficiently large, we have

$$ F_3 = \begin{pmatrix} F_{h,2} & F_{N,2} \\ F_{x,2}^* & F_z \end{pmatrix}, $$

where

$$ F_{h,2} = |S_R|^2 (t_r^2 |\mu|^2 \sin^2 \xi_r V_r V_r^* + R_h I)^\dagger, \quad (41) $$

$$ F_{x,2} = -|S_R|^2 t_r^2 \mu \cos \xi_r \sin \xi_r V_r q_r^* \dagger, $$

$$ V_r = \begin{pmatrix} \cos \theta_r q_r \\ \sin \theta_r q_r \end{pmatrix}, $$

$$ \mu = \tau_r + j. $$

$V_r$ is a $2N \times 1$ matrix. $F_{h,2}$ is a $2N \times 2N$ matrix. $F_{x,2}$ is a $2N \times N$ matrix. We can calculate the inverse of $F_3$ using eq. (22) because $V_r$ has the same property as $q_r$, that is, $V_r^* V_r = N$. Thus, eq. (10) is also valid for $V_r$, which leads to the inverse of $F_3$:

$$ F_3^{-1} = \frac{1}{R_h \lambda N |S_R|^2} \quad (42) $$

$$ \begin{pmatrix} \lambda N I - |\mu|^2 \sin^2 \xi_r V_r V_r^* & \sigma \mu \cos \xi_r \sin \xi_r V_r q_r^* \\ \sigma \mu^* \cos \xi_r \sin \xi_r q_r V_r^* & \sigma [\lambda N I - \sigma \cos^2 \xi_r q_r q_r^*] \end{pmatrix}^\dagger, $$

where

$$ \lambda = \sigma \cos^2 \xi_r + |\mu|^2 \sin^2 \xi_r + \frac{R_h}{N t_r^2}. \quad (43) $$
<table><caption>**Table 1.** First two derivatives of the steering matrices versus $\xi$. For $P_{3,a}$, $P_{3,h}$ and $P_{3,z}$ the expressions are obtained by replacing $q$ by $V$ for the first block in each matrix.</caption><thead><tr><th>P</th><th>E<sup>*</sup>(ξ)</th><th>∂E<sup>*</sup>(ξ) / ∂ξ</th><th>∂<sup>2</sup>E<sup>*</sup>(ξ) / ∂ξ<sup>2</sup></th></tr></thead><tbody><tr><td>P<sub>2,a</sub></td><td>(j sin ξ q<sup>*</sup> cos ξ q<sup>*</sup>)</td><td>E<sup>*</sup><sub>2,a</sub>(ξ + π/2)</td><td>-E<sup>*</sup><sub>2,a</sub>(ξ)</td></tr><tr><td>P<sub>2,h</sub></td><td>(j tan ξ q<sup>*</sup> q<sup>*</sup>)</td><td>j / cos<sup>2</sup> ξ (q<sup>*</sup> 0)</td><td>2 tan ξ ∂E<sup>*</sup><sub>2,h</sub>(ξ) / ∂ξ</td></tr><tr><td>P<sub>2,z</sub></td><td>(j q<sup>*</sup> cot ξ q<sup>*</sup>)</td><td>-1 / sin<sup>2</sup> ξ (0 q<sup>*</sup>)</td><td>-2 cot ξ ∂E<sup>*</sup><sub>2,z</sub>(ξ) / ∂ξ</td></tr></tbody></table>

requires $3N \times (2 \times 3N + 1) + 1$ operations (multiplications or additions, the cost of the unique division is counted as twice the cost of a multiplication, replace $3N$ by $2N$ for a projected cross-spectrum matrix), the product of $F_3^{-1}$ with $E_{3,a}$ being the most demanding part. The $3N$ values of this matrix can be saved while computing $P_{3,a}$, to be re-used for the first derivative. $3N$ operations (only one sine and one cosine for the whole matrix, whose values are already available) are necessary to initialize $E_{3,a}(\xi + \pi/2)$, $q$ being already available for the current wavenumber. There are finally $2 \times 3N - 1$ operations to multiply the constructed steering matrix with the saved matrix. Hence, the total additional cost of the first derivative is $3 \times 3N + 2$ operations (3.3 per cent of the cost of $P_{3,a}$ for 15 sensors, less for a higher number of sensors). The second derivative is more expensive because $P_{3,a}$ must be re-evaluated at $\xi + \pi/2$. On the whole, the cost of the function value and its two derivatives is slightly larger than twice the cost of $P_{3,a}$. The average number of equivalent $P_{3,a}$ evaluations observed for a typical run is 6.79 ± 2.91 per optimization (15 3C sensors recording during 1h58', processed between 0.5 and 10 Hz during 2'55" with 4 threads on a laptop - Intel Core i7-8650U). The same results are achieved with a nested grid search technique (initial step of 9°) by the means of 33.25 ± 11.13 evaluations per optimization (11'32" on the same machine and under the same conditions). Assuming that there is a common fixed CPU time (input/output, pre-processing, Fourier transforms,...) and a CPU time proportional to the number of $P$ evaluations, 76 and 94 per cent of the time is spent in computing $P$ for the Newton iterations and the nested grid search, respectively. The details of operation counting and CPU time are provided in Section 3 of Supporting Information.

## 3.6 Ellipticity correction

In the previous sections, we saw that the ellipticity estimators obtained by maximizing the beam power are biased by the incoherent noise level and its partition between the vertical and the horizontal components. A 3C processing provides only two independent equations (eqs 33 and 34) between the observations ($\hat{\xi}_h$, $\hat{\xi}_z$) and the unknowns ($\xi_r$, $R_h$ and $R_z$). An additional constraint to this problem can be found by considering the second derivative of the beam power in the radial direction at $\bar{k} = \bar{k}_r$. From eq. (46), we have

$$ \frac{\partial^2 P_{3,a}}{\partial k^2}(\bar{k}_r, \hat{\xi}_a) = c_a(\hat{\xi}_a, \xi_r)P_{3,a}(\bar{k}_r, \hat{\xi}_a) \left[ \frac{\partial^2 |B(\bar{k} - \bar{k}_r)|^2}{\partial k^2} \right]_{\bar{k}=\bar{k}_r}^\dagger, (49) $$

where

$$ \begin{aligned} c_a(\hat{\xi}_a, \xi_r) &= \frac{\rho(\hat{\xi}_a, \xi_r)}{\sin^2 \hat{\xi}_a + \sigma \cos^2 \hat{\xi}_a - \rho(\hat{\xi}_a, \xi_r)}, \\ \rho(\xi, \xi_r) &= \frac{(\sin \xi_r \sin \xi + \sigma \cos \xi_r \cos \xi)^2}{\lambda}, \end{aligned} \quad (50) $$

and $k$ is the absolute value of vector $\bar{k}$. Exactly the same relations are obtained for the derivative of $P_{2,a}$. In eq. (49), all quantities except

the concavity factor $c_a$ can be experimentally calculated, hence an experimental value for $c_a$ can be deduced. $\xi_r$ can be adjusted so that $c_a$ fits its experimental value ($\bar{c}_a$). The determination of the second derivative of $P_{3,a}$ and $|B(\bar{k})|^2$ versus $k$ is detailed in Appendix B. The second derivative of the beam pattern at $\bar{k} = \bar{k}_r$ is always negative, $P_{3,a}$ is always positive, $\hat{c}_a$ is thus positive because $P_{3,a}$ has a maximum at $\bar{k}_r$, namely, $\hat{\xi}_a$. The first two equations that link the three unknowns $\xi_r$, $R_h$ and $R_z$ are eqs (33) and (34) that can be transformed into:

$$ \begin{aligned} \frac{R_z}{N} &= \cos \xi_r \sin \xi_r \left( \frac{1}{\hat{e}_h} - \frac{1}{e_r} \right), \\ \frac{R_h}{N} &= \cos \xi_r \sin \xi_r (\hat{e}_z - e_r), \end{aligned} \quad (51) $$

where $e_r$, $\hat{e}_h$ and $\hat{e}_z$ are the tangents of $\xi_r$, $\hat{\xi}_h$ and $\hat{\xi}_z$, respectively. Substituting $\sigma$, $\rho$, $\frac{R_h}{N}$ and $\frac{R_z}{N}$ by their respective values in eq. (50), we have

$$ \hat{c}_a = c_a(e_r) = \frac{[(\hat{e}_a - \hat{e}_h)e_r + \hat{e}_h(\hat{e}_z - \hat{e}_a)]^2}{[(\hat{e}_a - \hat{e}_h)\hat{e}_a + \hat{e}_h(\hat{e}_z - \hat{e}_a)](\hat{e}_z - e_r)(e_r - \hat{e}_h)}, \quad (52) $$

which is a quadratic equation with unknown variable $e_r$. $\hat{c}_a$ is the tangent of $\hat{\xi}_a$. Therefore, there are two solutions to the problem without any way to discriminate: one corresponds to the true solution and the other one is an artefact. From the inequalities in eq. (38), function $c_a(e_r)$ is positive between $\hat{e}_h$ and $\hat{e}_z$. It has a single minimum at $\frac{\hat{e}_h + \hat{e}_z}{\hat{e}_a + \hat{e}_h}$ between two poles at $\hat{e}_h$ and $\hat{e}_z$. The value at the minimum is $\frac{4\hat{e}_h\hat{e}_a}{(\hat{e}_a - \hat{e}_h)(\hat{e}_a + \hat{e}_h(\hat{e}_z - \hat{e}_a))}$. Similar equations can be derived with the second derivatives of $P_h$ and $P_z$ in the radial direction at $\bar{k} = \bar{k}_r$. Replacing $\hat{c}_a$ by $\hat{e}_h$ or $\hat{e}_z$, the equation simplifies to

$$ \hat{c}_h = c_h(e_r) = \frac{\hat{e}_h(\hat{e}_z - \hat{e}_h)}{(e_r - \hat{e}_h)(\hat{e}_z - e_r)} \quad \text{and} \quad (53) $$

$$ \hat{c}_z = c_z(e_r) = \frac{e_r^2(\hat{e}_z - \hat{e}_h)}{\hat{e}_z(e_r - \hat{e}_h)(\hat{e}_z - e_r)}. \quad (54) $$

These two functions have the same properties as $c_a(e_r)$ with minima at $\frac{\hat{e}_h + \hat{e}_z}{2}$ and $\frac{2\hat{e}_h\hat{e}_z}{\hat{e}_h + \hat{e}_z}$, respectively. Their minimum values are both at $\frac{4\hat{e}_h}{\hat{e}_z - \hat{e}_h}$. They are also quadratic which bring two other pairs of solutions. $\hat{c}_h$, $\hat{c}_z$ and $\hat{c}_a$ must be greater than the minimum value otherwise the equations have no formal solutions. It might be tempting to combine eqs (53) and (54) to remove the ambiguity and to obtain an explicit expressions for $e_r$:

$$ \hat{e}_{h z} = \pm \sqrt{\frac{\hat{c}_z}{\hat{c}_h} \hat{e}_h \hat{e}_z}. \quad (55) $$

There is no sign ambiguity for the square root because $e_r$ must have the same sign as $\hat{e}_h$ and $\hat{e}_z$ (eq. 38). Two alternatives are also available by combining eqs (52) and (53) (or eq. 54) in the same way as in eq. (55) to calculate $\hat{e}_{h a}$ or $\hat{e}_{z a}$. Their analytical formulations are too complex to bring valuable information but they can be easily computed numerically. For RTBF method, the estimation of the
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Robustness of three ellipticity corrections in the case of randomly biased experimental quantities $\hat{e}_a, \hat{e}_h, \hat{e}_z, \hat{c}_a, \hat{c}_h$ and $\hat{c}_z$. (a) Ellipticity correction based on $\hat{c}_h$ and $\hat{c}_z$ ($\hat{e}_{hz}$, eq. 55). (b) Ellipticity correction based on $\hat{c}_h$ and $\hat{c}_a$ ($\hat{e}_{ha}$, combination of eqs 52 and 53). (c) Ellipticity correction based on $\hat{c}_z$ and $\hat{c}_a$ ($\hat{e}_{za}$, combination of eqs 52 and 54). The incoherent noise ratios $\frac{R_h}{N}$ and $\frac{R_z}{N}$ randomly vary from 0.1 to 10. The same thousand realizations are displayed in each panel with a drawing opacity of 10 per cent. The darkness is thus an indication of the density of dots.</figcaption></figure>

ellipticity ($\hat{\xi}_s$) is the square root of the product of $\hat{e}_h$ and $\hat{e}_z$ (eq. 37). By considering the second derivative of the beam power, eq. (55) shows that the incoherent noise influence can be corrected by a simple factor made of the ratio of the concavity factors observed for $P_h$ and $P_z$. For a single plane wave, $\hat{c}_h$ and $\hat{c}_z$ are always positive and they are unbiased, leading to a unique solution of the system of equations and a perfect removal of the incoherent noise influence. For ambient vibration wavefields where simultaneous wave arrivals may bring unwanted contributions to the beam power and its second derivative, $\hat{c}_h$, $\hat{c}_z$ and $\hat{c}_a$ can be biased or even negative. Therefore, the system of three equations has no formal solution despite that eq. (55) always brings a real solution if $\hat{c}_h$ and $\hat{c}_z$ have the same sign.

Fig. 5 compares the robustness of the three ellipticity corrections ($\hat{e}_{hz}, \hat{e}_{ha}$ and $\hat{e}_{za}$) in the case of imperfect experimental quantities $\hat{e}_a, \hat{e}_h, \hat{e}_z, \hat{c}_a, \hat{c}_h$ and $\hat{c}_z$ for random $\frac{R_h}{N}$ and $\frac{R_z}{N}$ varying from 0.1 to 10. $\hat{e}_h$ and $\hat{e}_z$ are modified randomly with a maximum error of $\pm 50$ per cent compared to their theoretical values. $\hat{e}_a$ is computed with eq. (36) which remains valid even for noisy experimental data. The random perturbations of $\hat{c}_a, \hat{c}_h$ and $\hat{c}_z$ are limited to $\pm 10$ per cent. All perturbations are generated with a random uniform generator over a logarithmic scale. Thousand realizations are plotted in Fig. 5 for true ellipticities between 1 and 89° (step of 1°). $\hat{e}_{hz}$ offers the best correction with a remarkable stability. For the chosen random intervals, a maximum of 12° deviation is observed at 45°, vanishing at both ends of the range. The angular scale ($\xi_r$) in Fig. 5 suggests a perfect robustness at 0° and $\pm 90$° whereas the relative error on the ellipticity ($e_r$) is approximately independent of the true ellipticity and around 50 per cent (not shown). It is an effect of the angular transformation. The other two corrections ($\hat{e}_{ha}$ and $\hat{e}_{za}$) are not as efficient with unstable behaviours for some realizations, in particular

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Robustness of ellipticity corrections based on the minimization of $\Phi$ functions in the case of randomly biased experimental quantities $\hat{e}_a, \hat{e}_h, \hat{e}_z, \hat{c}_a, \hat{c}_h$ and $\hat{c}_z$: (a) $\Phi_{hz}$, (b) $\Phi_{hza}$, (c) $\Phi_{ha}$ and (d) $\Phi_{za}$. The same thousand realizations as in Fig. 5 are displayed with the same options. $\phi_m$ is equal to 100.</figcaption></figure>

for true ellipticities lower than 20° and greater than 70° for $\hat{e}_{ha}$ and $\hat{e}_{za}$, respectively.

According to Fig. 2(a), $\hat{\xi}_h$ is almost unbiased around 0° and the beam power is higher for low than for high absolute ellipticity values. In experimental conditions, low beam power peaks have more chance to be distorted by neighbouring peaks than stronger ones. The distortion might be such that at $\hat{k} = \hat{k}_r$ the concavity is no longer negative. Hence, $\hat{c}_h$ is likely to be biased or even negative at high absolute ellipticity values. The exact opposite situation is encountered for $\hat{\xi}_z$ and $\hat{c}_z$. Therefore, $\hat{e}_{hz}$ cannot be used at low and high absolute ellipticity values. The help of $\hat{e}_{ha}$ or $\hat{e}_{za}$ is thus required precisely where these two estimators are not robust (Figs 5b and c). Alternatively, we could solve the ambiguity of eqs (52), (53) and (54) in the least-square sense by searching for the absolute minimum of either of the positive functions:

$$ 
\begin{align} 
\Phi_{hz}(e_r) &= \left( \frac{c_h(e_r)}{\hat{c}_h} - 1 \right)^2 + \left( \frac{c_z(e_r)}{\hat{c}_z} - 1 \right)^2, \notag \\ 
\Phi_{hza}(e_r) &= \left( \frac{c_h(e_r)}{\hat{c}_h} - 1 \right)^2 + \left( \frac{c_z(e_r)}{\hat{c}_z} - 1 \right)^2 + \left( \frac{c_a(e_r)}{\hat{c}_a} - 1 \right)^2, \notag \\ 
\Phi_{ha}(e_r) &= \left( \frac{c_h(e_r)}{\hat{c}_h} - 1 \right)^2 + \left( \frac{c_a(e_r)}{\hat{c}_a} - 1 \right)^2 \text{ or} \notag \\ 
\Phi_{za}(e_r) &= \left( \frac{c_z(e_r)}{\hat{c}_z} - 1 \right)^2 + \left( \frac{c_a(e_r)}{\hat{c}_a} - 1 \right)^2. \tag{56}
\end{align} 
$$

These four functions can be reduced to fractions of polynomials: degree 4 over 4. Their first derivatives are also fractions of polynomials: degree 6 over 8. There is no pole between $\hat{e}_h$ and $\hat{e}_z$. The four functions are thus continuous on the interval of interest. According to the degree of the numerator of their first derivatives, there are at most 3 minima (corresponding to 5 real roots of the derivative). The roots of the first derivative can be solved by constructing the companion matrix of the numerator. Its eigenvalues are all the complex roots of the polynomial. Roots with a non null imaginary part are discarded. They are rounded to 0.0001° (after conversion to angular
ellipticity) to avoid multiplicity within numerical precision. Roots outside the interval $\hat{e}_h$ to $\hat{e}_z$ are ignored. If there are more than one minimum, only the lowest is kept. To avoid ambiguous cases, there must be a minimum factor ($\phi_m$) between the value of the function at the best minimum and at all other minima.

In Fig. 6, these four functions are analysed with the same random realizations as in Fig. 5 and with $\phi_m$ equal to 100. The correction is effective on most cases. The number of outliers is clearly reduced compared to $\hat{e}_{ha}$ and $\hat{e}_{za}$. $\Phi_{hz}$ is slightly less stable than $\Phi_{hz}$. From a set of 1000 random realizations, it is possible to measure various indicators: the maximum error, the number of accepted values (some realizations are rejected because of $\phi_m$ condition or if $|e_h| > |e_z|$) and the number of values shifted by less than 5° and by more than 10°. To get robust values obtained after a random process, the operations are repeated 20 times. A mean and a standard deviation is calculated for each of these indicators. They are summarized in Table 2 for all ellipticity estimators and for several values of $\phi_m$. $\hat{e}_{hz}$ produces the lowest number of outliers but at the price of providing a good correction (<10° error) for only 72 per cent of the cases, although $\Phi_{hz}$ has slightly more erroneous corrections but good corrections are available for 95 per cent of the cases. Though the maximum error is kept within a small interval (<15° error). Increasing $\phi_m$ from 2 to 1000 reduces the number of accepted corrections by approximately 15 per cent for $\Phi_{hz}$ (40 per cent for $\Phi_{ha}$ and $\Phi_{za}$) but the number of bad corrections (>10° error) is almost not affected except for $\Phi_{ha}$ and $\Phi_{za}$ (reduction by a factor of 2).

From the considerations hereinabove, it is suggested to correct the ellipticity with $\Phi_{hz}$ if $\hat{e}_h$ and $\hat{e}_z$ are positive. $\phi_m$ equal to two in case of multiple minima is retained because higher values are not sufficiently effective to discard erroneous cases. If either $\hat{e}_h$ or $\hat{e}_z$ is positive but not both, $\Phi_{ha}$ or $\Phi_{za}$ is used with $\phi_m$ equal to 100. Although $\hat{e}_{hz}$ is faster to compute than searching for the minimum of $\Phi_{hz}$, it is not chosen because it is less tolerant of bias in experimental quantities.

In the next sections, synthetic and experimental examples are provided to illustrate the effects of incoherent noise on ellipticity results and their corrections with the proposed algorithm.

## 4 A SYNTHETIC AMBIENT VIBRATION EXAMPLE

The same synthetic data set as in Wathelet *et al.* (2018) is considered which provides signals during 6 hr at 38 receivers for model M2.1 (SESAME European Project, 'Site EffectS assessment using AMbient Excitations', Bard *et al.* 2004). It is made of a soft layer over a rigid half-space with a strong velocity contrast. Signals are available at a sampling rate of 200 Hz. However the sources are providing excitation of the structure up to 50 Hz. Only the results obtained with subset *array 2* are presented hereinafter (14 sensors, 75 m aperture). The geometry of this array shown in Fig. S2 is the same as in Poggi & Fäh (2010, fig. 7).

In Wathelet *et al.* (2018), sets of 4N blocks were used to compute the cross-spectral matrices of size $2N \times 2N$. Experience has shown that such block sets are not large enough to get a stable inverse for cross-spectral matrices of size $3N \times 3N$. Therefore, processing sets of at least 6N blocks for ARDS is recommended. For a fair comparison between methods, the same number of blocks is chosen everywhere, even if smaller numbers are still acceptable for VBF, RTBF and ARTBF. The length of the blocks is fixed to 50 periods instead of 100 in Wathelet *et al.* (2018). This reduction does not alter significantly the results and reduces the minimum required

recording time. For each frequency and block set, a maximum of N local maxima of the beam power, for which the amplitude is greater than 1 per cent of the global maximum, is searched over the $k_x$, $k_y$ and $\xi$ domain using the technique described in Section 3.5. The other details of the process are exactly the same as in Wathelet *et al.* (2018).

In the following discussion, the modes of Rayleigh waves are numbered $r_0, r_1, r_2, ...$ for the fundamental mode and the higher order modes. For Love waves, the same convention is adopted: $l_0, l_1, l_2, ...$.

In Figs 7(a)–(d), one histogram per frequency is built from the velocity values produced by methods VBF, RTBF, ARTBF and ARDS, respectively. Mode $r_0$ is correctly retrieved for all four methods above 2.3 Hz, even if denser aliasing patterns are observed for increasing wavenumbers, in particular above $k_{max}/2$ (Wathelet *et al.* 2008). RTBF and ARTBF do not provide a sharper image of mode $r_0$ than a single component processing like VBF, whilst ARDS delivers a more precise and accurate dispersion curve between 2.1 and 6 Hz. Additionally, it deviates from the true dispersion curve at a slightly lower frequency than VBF (2.3 Hz) and RTBF (2.2 Hz). ARTBF shows artefacts between 2 and 4 Hz that blur the dispersion curve below 2.5 Hz. They can be attributed to mode $l_0$. Under the same conditions, RTBF and ARDS are not influenced at all by the presence of Love waves, that is, the separation between Love and Rayleigh waves on the horizontal components is effective. Below a limit that varies from 2.1 to 2.3 Hz, the velocity is clearly overestimated everywhere, as expected for wavelengths greater than four or five times the array aperture.

Mode $r_1$ is much better retrieved for all 3C methods than for VBF (between 4 and 8 Hz). Like for mode $r_0$, there is a slightly higher precision for ARDS than for other methods. Modes $r_1$ and $r_2$ look like artificially connected between 7.3 and 8.3 Hz with no specific advantage for one method. Other modes are also visible but partially blurred by aliasing patterns, probably preventing any reliable picking without the help of the theoretical curves shown in Fig. 7. Aliasing patterns are notably reduced when considering, in their order, Figs 7(a)–(d). ARDS has the cleanest image with no spurious velocities below 8 Hz between modes $r_0$ and $r_1$, in particular at 3.75 and 6 Hz.

The raw ellipticity curves are observed by constructing one histogram at each frequency for all $\xi$ values that maximize the beam power. They are obtained jointly with the dispersion curves for RTBF, ARTBF and ARDS (Figs 8a, b and c, respectively). One ellipticity value corresponds to one velocity value and vice versa. There is no case where the ellipticity cannot be determined. For a clearer view of the ellipticity curves of each mode, the dispersion curves are manually picked where they are obviously visible and unambiguously separable from other modes. The velocity values located at less than 20 per cent from the picked curves are selected to build the histograms shown in Figs S3 and S4 for ARTBF and ARDS, respectively.

The most striking improvement brought by the all-component ellipticity steering is illustrated by the difference between RTBF and the other two methods for ellipticity values close to zero. It is impossible to get values between -15° and +15° with RTBF: the results are shifted towards higher absolute values. In eq. (35), for the neighbourhood of zero ($|\xi| < 15°$), the squared sine of the numerator does not exceed 0.067. The squared cosine of the denominator is close to 1 and at least 14 times greater than the squared sine. Therefore, even a very small amount of incoherent noise on the horizontal components can shift the observations. A similar behaviour is noted for angular ellipticities close to $\pm 90°$.
<table><caption>**Table 2.** Comparison of the ellipticity corrections for random perturbations of $\hat{e}_a$, $\hat{e}_h$, $\hat{e}_z$, $\hat{c}_a$, $\hat{c}_h$ and $\hat{c}_z$. The perturbations are the same as in Fig. 5. The incoherent noise ratios $R_h/N$ and $R_z/N$ randomly vary from 0.1 to 10. 1000 random realizations are calculated for true ellipticities between 1° and 89° (step of 1°, 89 000 ellipticity corrections). $\phi_m$ is the minimum factor between the value of the function at the best minimum and at all other minima. For each estimator, the following indicators are available: the maximum error, the number of accepted values and those falling below a 5° above a 10° deviation. For the sake of robustness, the process is repeated 20 times. Means and standard deviations (shown after $\pm$) are computed over 20 samples for each case.</caption><thead><tr><th>Estimator</th><th>φ<sub>m</sub></th><th>Max. error</th><th>Total accepted</th><th>Error &lt;5°</th><th>Error &gt;10°</th></tr></thead><tbody><tr><td>ε̂<sub>hz</sub></td><td>-</td><td>12 ± 0</td><td>64556 ± 1023</td><td>54036 ± 1099</td><td>460 ± 85</td></tr><tr><td>ε̂<sub>ha</sub></td><td>-</td><td>87 ± 0</td><td>60740 ± 679</td><td>34965 ± 673</td><td>12676 ± 502</td></tr><tr><td>ε̂<sub>za</sub></td><td>-</td><td>88 ± 0</td><td>56162 ± 952</td><td>31902 ± 637</td><td>13562 ± 598</td></tr><tr><td>Φ<sub>hz</sub></td><td>2</td><td>15 ± 1</td><td>87588 ± 172</td><td>64086 ± 607</td><td>2985 ± 317</td></tr><tr><td></td><td>10</td><td>16 ± 1</td><td>82525 ± 568</td><td>61316 ± 971</td><td>2594 ± 248</td></tr><tr><td></td><td>100</td><td>16 ± 1</td><td>77276 ± 808</td><td>57057 ± 792</td><td>2455 ± 237</td></tr><tr><td></td><td>1000</td><td>15 ± 1</td><td>75354 ± 661</td><td>55351 ± 738</td><td>2626 ± 231</td></tr><tr><td>Φ<sub>hza</sub></td><td>1000</td><td>78 ± 10</td><td>71620 ± 642</td><td>51625 ± 714</td><td>3324 ± 226</td></tr><tr><td>Φ<sub>ha</sub></td><td>2</td><td>87 ± 1</td><td>84210 ± 260</td><td>54327 ± 479</td><td>9585 ± 426</td></tr><tr><td></td><td>10</td><td>87 ± 1</td><td>71670 ± 680</td><td>47457 ± 930</td><td>6874 ± 365</td></tr><tr><td></td><td>100</td><td>87 ± 1</td><td>57720 ± 879</td><td>38414 ± 797</td><td>5274 ± 318</td></tr><tr><td></td><td>1000</td><td>86 ± 3</td><td>49916 ± 822</td><td>32795 ± 655</td><td>4757 ± 272</td></tr><tr><td>Φ<sub>za</sub></td><td>2</td><td>86 ± 3</td><td>84317 ± 214</td><td>53806 ± 720</td><td>9520 ± 425</td></tr><tr><td></td><td>10</td><td>87 ± 1</td><td>71715 ± 719</td><td>46650 ± 693</td><td>6836 ± 275</td></tr><tr><td></td><td>100</td><td>87 ± 1</td><td>57296 ± 753</td><td>37616 ± 722</td><td>5167 ± 331</td></tr><tr><td></td><td>1000</td><td>87 ± 1</td><td>48372 ± 1231</td><td>31460 ± 1071</td><td>4505 ± 276</td></tr></tbody></table>

The gaps around zero and around $\pm 90°$ are mainly controlled by the horizontal and the vertical incoherent noise ratios, respectively. Such effects completely disappear when considering $\hat{\xi}_a$ (eq. 32) as demonstrated by the results in Figs 8(b) and (c) where the full range of ellipticity is covered without notable bias below 8 Hz (compared to the theoretical curves). Above 8 Hz, aliasing patterns are blurring the image: more than one mode of propagation is encountered which makes the situation more complex, that is, for a particular wave, the other waves with distinct wavenumbers and ellipticities can be partly viewed as incoherent noise. The perceived incoherent noise is thus increasing with the frequency which leads to a convergence of the ellipticities towards $\pm 45°$ for RTBF as already noted while commenting eqs (19) and (35).

ARTBF and ARDS provide a remarkable image of the ellipticity at least up to mode $r_5$ where the rapid transition between 14 and 21 Hz is even detected (Figs S3f and S4f). It can be verified that the higher density patches on the dispersion curves have their equivalent on the ellipticity curves and vice versa. Moreover, they approximately match with the same mode in both domains. Nonetheless, there are two main differences between ARTBF and ARDS. (1) Between 1.75 and 4 Hz, a lot of values close to $\pm 90°$ encountered for ARTBF do not appear on ARDS results. They are related to the velocities located in the neighbourhood of mode $l_0$ noted in Fig. 7(c) between 1.75 and 4 Hz. An absolute angular ellipticity of $90°$ is consistent with Love waves. Between 2.4 and 4 Hz, the high absolute ellipticities can be easily filtered out to clean the dispersion picture. However below 2.4 Hz, Rayleigh and Love waves have approximately the same ellipticity and their dispersion curves merge which prevents any unambiguous separation. (2) Below 5 Hz, $r_0$ ellipticity curve from ARDS is more precise than that from ARTBF, which allows a proper differentiation of modes $r_0$ and $r_1$ at 4 Hz. Mode $r_1$ is equally retrieved by both methods. This case is favourable because the dispersion curves of the two modes are sufficiently far apart (relative to the array resolution) when the two ellipticity curves are kissing each other.

Between modes $r_1$ and $r_2$, spurious velocity estimations are observed between 7.3 and 8.3 Hz as already mentioned hereinabove.

The ellipticity provides a direct explanation: modes $r_1$ and $r_2$ are kissing each other in terms of ellipticity around 7 Hz. The ellipticity difference between modes $r_0$ and $r_1$ at 4 Hz and between modes $r_1$ and $r_2$ at 8.3 Hz are $18°$ and $16°$, respectively. The wavenumber difference (between the corresponding modes) at 4 and 8.3 Hz are 0.056 and 0.075 rad m$^{-1}$, respectively. Below 8.3 Hz, $r_1$ and $r_2$ are closer down to $1.5°$ which makes the separation impossible. Above 8.3 Hz, the distance between modes (in velocity and ellipticity terms) is almost similar as the distance at 4 Hz: the modes are well identified by ARTBF and ARDS. The relative power between modes may also play a role in the separation capabilities.

For this synthetic data set, the ellipticity correction introduced in Section 3.6 does not modify significantly the ellipticity image as shown in Figs S5 and S6 by comparison with Figs 8(c) and S4, respectively. The fundamental mode has less values far from the theoretical curve around 20 Hz after correction, which remains a minor improvement. One interest of the proposed correction is also a separate assessment of the incoherent noise levels for horizontal and vertical components. In Fig. 9(a) and (b), experimentally calculated $R_h/N$ and $R_z/N$ are as low as 0.01 for frequencies below 4 Hz. They tend to hundred times this value at high frequencies. Aliasing at high frequencies explains the progressive increase of the noise with the frequency. Below 3.5 Hz, there is apparently an ambiguity in the determination of $R_h/N$ and $R_z/N$. Given the very low level of incoherent noise, it has a minor impact on the ellipticity correction: a difference of about $0.5°$ is observed for the two solutions. $\delta$ shown in Fig. 9(c) is almost flat and close to zero below 5 Hz except at 3.8 Hz where the true ellipticity vanishes for $r_0$. Above 5 Hz, the ellipticity uncertainty increases together with the average $\delta$. Considering mode $r_1$ alone in Fig. S6(b), the flat part extends up to 7 Hz except when the true ellipticity approaches $0°$ or $\pm 90°$. For both modes, the retrieved ellipticity ($\hat{\xi}_a$) is almost unbiased when $\delta$ is close to 0. $\delta$, being the deviating term in eq. (32), is thus a good indicator of the accuracy of $\hat{\xi}_a$. Fig. 10 shows only the dispersion and ellipticity results for which $|\delta|$ is less than 0.1. Among notable differences compared to previous figures is the proper discarding of values around the dispersion curve of $r_0$ and above 20 Hz, better highlighted in Fig. S7(a) where modes
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison of the dispersion curves obtained by processing the 6-hr-long signals simulated by Wathelet *et al.* (2018) on array 2 for the four beam powers: (a) VBF (*P*<sub>z</sub>), (b) original RTBF (*P*<sub>2,s</sub>), (c) all-component ellipticity steering RTBF (ARTBF, *P*<sub>2,a</sub>) and (d) ARDS (*P*<sub>3,a</sub>). The colour scale [batlow from Crameri *et al.* (2020)] represents the probability density of the velocity histograms built for each frequency. The black straight lines are the theoretical array limits defined by Wathelet *et al.* (2008), from left to right: *k*<sub>min</sub>/2 (plain line), *k*<sub>min</sub> (dotted line), *k*<sub>max</sub>/2 (dotted line) and *k*<sub>max</sub> (dashed line). The black dotted curves are the theoretical Rayleigh dispersion curves for the fundamental (without symbols) and the first five higher order modes (in order, crosses, circles, triangles, stars and diamonds). In panel (c), the blue dotted curve is the theoretical Love dispersion curve for the fundamental mode.</figcaption></figure>

are split apart. Even if the proposed correction is not perfect, this
example demonstrates that the additional noise indicators can help
in identifying reliable parts of the curves.

In the results discussed hereinabove, all horizontal components
are perfectly oriented and all sources are located far from the sen-
sors, therefore the incoherent noise levels can be very low as shown
in Fig. 9, leading to negligible deviations of ellipticity estimators.
To test the effectiveness of the ellipticity correction, random noise
is deliberately added to the simulated signals, either on the vertical
or on the horizontal components, or even on both. Only the hori-
zontal perturbation is presented hereinafter. After the generation of
42 (14 sensors × 3 components) random signals between −1 and
1, their spectra are whitened and multiplied by 5 × 10<sup>-5</sup> before
being added to the simulated signals. The spectrum of a typical
horizontal signal and of the added noise is shown in Fig. S8. The
signal spectrum is greater than the noise spectrum only between
2 and 3 Hz which represents a very noisy situation. The modified
signals are processed with ARDS in the same way as the previous
results. Fig. 11 shows the dispersion and the ellipticity curves with-

out and with the ellipticity correction. The amount of incoherent
noise reduces the achievable results to a small part of modes r₀ and
r₁ (2.2–6.5 Hz). The raw ellipticity is clearly biased towards ±90°
for mode r₀ in accordance with eq. (32). The correction works fine
below 4 Hz whilst being visible but less efficient up to 6 Hz: -33°
is found instead of -28°. Mode r₁ between 4 and 5 Hz is well
delineated even without a correction. Above 5 Hz, the correction is
necessary to obtain a valid ellipticity up to 6 Hz. Above 6 Hz, the
correction is no longer effective.

In Fig. 8, results from array processing are compared to RayDec,
a single station processing proposed by Hobiger *et al.* (2009). The
mean and the standard deviation are computed from the curve of
each station considered individually. As recommended, ten cycles
are stacked after filtering with a relative bandwidth of 0.2. The ob-
tained unsigned curve matches with the true ellipticity between 2.4
and 3.3 Hz. At the ellipticity peak, the smooth RayDec curve under-
estimates the true ellipticity in accordance with Hobiger *et al.* (2009)
whose best results deviated from theoretical ellipticities above 80°.
Above 3.3 Hz, RayDec does not support a vanishing horizontal
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of the ellipticity curves obtained simultaneously with the dispersion curves shown in Fig. 7: (a) RTBF ($P_{2,s}$), (b) ARTBF ($P_{2,a}$) and (c) ARDS ($P_{3,a}$). The light blue curves are the mean (plain line) and one standard deviation (dashed lines) obtained for RayDec method over all stations of array 2 (Hobiger et al. 2009). They are shown twice for positive and negative ellipticities because RayDec methods provide only the absolute value. The black dotted curves are the theoretical Rayleigh wave ellipticity curves for the first six modes with the same legend as in Fig. 7.</figcaption></figure>

component probably being influenced by higher order modes or other waves (body and Love waves). Considering the negative duplicated curve, the absolute value of the fundamental mode is correctly retrieved above 8 Hz. Between 4 and 6 Hz, mode $r_1$ has probably an influence. Array results provide signed values and they slightly

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Incoherent noise indicators obtained after calculating the ellipticity corrections for ARDS results in Fig. 8(c): (a) incoherent noise ratio for vertical components $R_\xi/N$, (b) incoherent noise ratio for horizontal components $R_b/N$ and (c) difference between $R_\xi/N$ and $R_b/N$ ($\delta$).</figcaption></figure>

better follow the true ellipticity curve around 2 Hz. Therefore, single station results can be used to validate the array results only for moderate ellipticity values (from 40° to 75° in this case) but not in the peak vicinity. Single station results are also influenced by the amount of incoherent noise on the horizontal components but not in the same proportion as uncorrected array ellipticities (Fig. 11). The corrected array ellipticity is more accurate and it is available over wider frequency range than RayDec ellipticity.

This synthetic example attests that ellipticity curves are better retrieved by ARTBF or ARDS than with RTBF, in particular around extremal values (0° and ±90°). Though ARTBF may not properly separate Rayleigh and Love waves like RTBF and ARDS. ARTBF and ARDS are more exposed to a confusion between the two polarizations than RTBF because they support the vanishing of the vertical component. Apparently, the inversion of the raw 3C cross-spectrum matrix offers an advantage to ARDS for the rejection of Love contributions. The proposed ellipticity correction proves to be efficient in a noisy context even if it is not always optimal. Nevertheless the indirect indicators remain valuable tools for the appraisal of ellipticity curves. In the next sections, several experimental data sets are analysed with the same methodology.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Histograms of dispersion and ellipticity curves for ARDS results with corrected ellipticity and after selecting only the values for which $|δ|$ (the absolute difference between the vertical and the horizontal incoherent noise ratios) is less than 0.1: (a) dispersion curves and (b) ellipticity curves. The black dotted curves have the same meaning as in Figs 7 and 8. The colour scales are the same as in Figs 7 and 8.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Histograms of dispersion and ellipticity curves for ARDS results with added incoherent noise on the horizontal components: (a) dispersion, (b) raw ellipticity and (c) corrected ellipticity curves. The blue and black dotted curves have the same meaning as in Figs 7 and 8. The colour scales are the same as in Figs 7 and 8.</figcaption></figure>

# 5 EXPERIMENTAL CASES

*InterPACIFIC* project (Garofalo et al. 2016a, b) investigated three sites (Mirandola in Italy, Grenoble and Cadarache in France) with ambient vibration arrays and borehole techniques. The publicly available data sets are used hereinafter as experimental test cases for the proposed beamforming methods. All ambient vibration measurements are sharing the same instrumentation (Güralp broadband CMG-6TD with integrated digitizers) and the same array geometries: a central station and two circles of 7 stations with a ratio of 3 between radii. The theoretical array response of this configuration is commented in Section 4 of the supplementary materials. The aliasing limit is $\bar{k}_{max}$ instead of $k_{max}$ specified in Wathelet et al. (2008). The recordings last from 1 to 2 hr. The distribution of array apertures varies among sites in relationship with the size of the shallow structure. These deployments are described with all details in Garofalo et al. (2016a). Data are processed with VBF, RTBF, ARTBF and ARDS with exactly the same parameters as for the synthetic case.

a regular shape from 0.8 to 4 Hz for ARDS (Fig. 12d). Mode $r_0$ is consistent with the fundamental mode reported in Garofalo et al. (2016a, fig. 3; see Fig. S11 for a formal comparison). RTBF and ARTBF results are pretty similar to those of ARDS except a discontinuous and more dispersed mode $r_1$ below 1.75 Hz for RTBF and a patch between 0.8 and 1 Hz just below $r_0$ for ARTBF. Like for the synthetic case, a contamination by Love waves can be suspected. However, the comparison with the Love dispersion curve obtained with a transverse high resolution beamforming (Fäh et al. 2008, Fig. S12 for details) does not firmly support this hypothesis. VBF overestimates the velocity below 1.4 Hz, probably because the separation between the two modes $r_0$ and $r_1$ is not as effective as for a 3C processing. Unlike all 3C methods, VBF is not able to catch the higher mode except between 2.5 and 4 Hz.

RTBF and uncorrected ARTBF and ARDS ellipticities are compared with the average RayDec ellipticity curve in Fig. 13. RTBF and RayDec fairly agree between 0.7 and 1.5 Hz as already observed by Wathelet et al. (2018). Uncorrected ARTBF and ARDS ellipticities are similar and they have a more realistic and continuous

## 5.1 Mirandola, Italy

It is a flat and relatively quiet area located on top of alluvial deposits (sand, silts and clays). At a depth between 50 and 150 m, stiffer layers are encountered (Garofalo et al. 2016a). Data sets labelled *PC1* to *PC4* (apertures from 30 to 810 m) are processed. Array *PC4* which has an aperture of 810 m is analysed in details hereinafter. The results of the other data sets are available in Fig. S10. The velocity histograms for the four methods applied to array *PC4* are shown in Fig. 12. Modes $r_0$ and $r_1$ are well individualized with
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Comparison of the dispersion curves retrieved for Mirandola array PC4 (aperture = 810 m) with (a) VBF, (b) RTBF, (c) ARTBF and (d) ARDS. The black straight lines are the theoretical array limits, from left to right: $k_{min}/2$ (plain line), $k_{min}$ (dotted line), $k_{max}/2$ (dotted line) and $\bar{k}_{max}$ (dashed line). The purple curves found in all panels are modes manually picked from panel (d): plain line for $r_0$ and dashed line for $r_1$. The blue curves are picked for the smaller array PC3 (Fig. S10). The dashed black curve in panel (c) is manually picked from a high resolution transverse beam forming of the same data, supposed to be mode $l_0$ (Fig. S12). The colour scales are the same as in Fig. 7.</figcaption></figure>

shape between 0.7 and 2.4 Hz, with a remarkable transition from prograde to retrograde. However, the comparison with RayDec is not as good as for RTBF below 1.5 Hz (about 6° error at 1 Hz for ARDS), ARTBF being slightly further than ARDS by approximately 4° at 1 Hz. Between 0.8 and 1 Hz, ARTBF results are likely to be influenced by Love waves which have an ellipticity of 90°. The shift direction is consistent with a Love wave influence as observed for the synthetic case in Fig. 8(b). Mode $r_1$ is clearly delineated even for the transition around ±90°, unlike the results produced by RTBF. The comparison with the results of the smaller array PC3 is rather good between 1.7 and 2.2 Hz for $r_0$ and $r_1$. Above 2.2 Hz the results of PC4 are probably distorted towards zero for both modes.

Fig. 14 shows the corrected ellipticity and the incoherent noise indicators for ARDS. Below 1.3 Hz, the correction decreases the ellipticity by approximately 5°, making the array results closer to RayDec results (1° difference at 1 Hz). Nevertheless, below 1 Hz, corrected ARDS and RayDec ellipticities diverge, the former approaching 90° while the latter reaches a maximum of 68° at 0.75 Hz. Given the similarities with the synthetic case when the ellipticity of the fundamental mode is close to ±90°, the ellipticity deduced from the proposed array processing is probably more reliable than the ellipticity obtained after a single station processing. The incoherent noise on the vertical components is much lower than on the

horizontal components below 1.3 Hz which explains the difference between corrected and uncorrected ellipticities. The corrected ellipticity between 2.3 and 3.1 Hz is significantly lowered compared to the uncorrected ellipticity and is in agreement with the ellipticity found for the smaller array PC3 (Fig. S10a). Above 3 Hz, the difference ($\delta$) between $R_c/N$ and $R_b/N$ is positive and monotonously increasing with the frequency: a concentration around zero is found for the ellipticity. The maximum frequency for which reliable ellipticity values can be retrieved is thus around 3 Hz which corresponds to approximately two times $\bar{k}_{max}$, even if the dispersion curve can be delineated for higher frequencies.

## 5.2 Grenoble, France

The site is in the middle of an industrialized area, surrounded by two highways and two rivers. The first tens of meters are characterized by recent stiff alluvial materials. The bedrock is encountered at depths from 500 to 800 m (Garofalo et al. 2016a). Data sets labelled PC1, PC2, PC3 and PC5 (apertures from 30 m to 810 m) are processed. Array PC3 which has an aperture of 270 m is analysed in details hereinafter. The results of the other data sets are available as supplementary materials in Fig. S13. The velocity histograms for the four methods applied to array PC3 are shown in Fig. 15. The same dispersion curve is retrieved for all methods from 1 to 10 Hz,
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Comparison of the ellipticity curves obtained simultaneously with the dispersion curves shown in Fig. 12: (a) RTBF, (b) ARTBF and (c) ARDS. The light blue curves are the mean (plain line) and one standard deviation (dashed lines) obtained for RayDec method over all stations of Mirandola array *PC4* (Hobiger *et al.* 2009). They are shown twice for positive and negative ellipticities because RayDec methods provide only the absolute value. The purple curves found in all panels are modes manually picked from panel (c): plain line for $r_0$ and dashed line for $r_1$. The dark blue curves are picked for the smaller array *PC3* (Fig. S10). The colour scales are the same as in Fig. 8.</figcaption></figure>

except for ARDS where the velocity is overestimated between 1.1 and 2 Hz compared to all other methods. A higher mode can be identified for RTBF (between 1.2 and 2 Hz) and ARTBF (between 0.85 and 1.7 Hz). On the contrary, ARDS does not detect this higher order mode correctly which induces a distortion of mode $r_0$ between 1.1 and 1.8 Hz where the distance between the two modes is minimum and close to the array limit $k_{min}$. For the larger array (PC5), the higher order mode is correctly retrieved and it is

<figure><img src="image_15.png" /><figcaption>**Figure 14.** ARDS results with corrected ellipticity for array PC4 (Mirandola): (a) ellipticity curves, (b) incoherent noise ratio for vertical components $R_\xi/N$, (c) incoherent noise ratio for horizontal components $R_\eta/N$ and (d) difference between $R_\xi/N$ and $R_\eta/N$ (δ). The light blue curves are the same as in Fig. 13. The purple curves are picked on the uncorrected results (Fig. 13) while the black curves are related to the corrected ellipticity. The dark blue curves are picked for the smaller array PC3 (corrected ellipticity, Fig. S10). The colour scales are the same as in Figs 8 and 9.</figcaption></figure>

compatible with the observations made for array PC3 (Fig. S13a). Despite this detail, the ARDS curve fits quite well the average curve found by (Garofalo *et al.* 2016a, fig. 9). A comparison can be found in Fig. S14.

The ellipticity obtained with RTBF is in agreement with RayDec above 0.9 Hz (Fig. 16). Below 1.4 Hz, ARTBF uncorrected ellipticities are compatible with the ellipticities obtained with the larger array (PC5) and they are diverging from RayDec estimation below 1.5 Hz. ARDS values are greater than those from ARTBF, leading to a difference of 25° between RayDec and ARDS at 1 Hz. There is no agreement between arrays PC5 and PC3 for ARDS below 1.8 Hz, which is precisely the range where the velocity is also dubious. ARDS results in this range can be reasonably considered as inaccurate. In Fig. S13(a), the results obtained for array PC5 also deviate from RayDec curve below 1.4 Hz. The array ellipticity delineates a tendency towards a prograde-retrograde transition which is a critical zone for RayDec as shown for the previous cases.
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Comparison of the dispersion curves retrieved for Grenoble array PC3 (aperture = 270 m) with (a) VBF, (b) RTBF, (c) ARTBF and (d) ARDS. The black straight lines are the theoretical array limits, from left to right: $k_{min}/2$ (plain line), $k_{min}$ (dotted line), $k_{max}/2$ (dotted line) and $\bar{k}_{max}$ (dashed line). The purple curve found in all panels is manually picked from panel (d). The red and the blue curves are picked for the larger (PC5) and the smaller (PC2) arrays, respectively (Fig. S13). The colour scales are the same as in Fig. 7.</figcaption></figure>

At 1 Hz, a significant deviation is observed while being in the optimum wavenumber range from the array response point of view (Fig. S13a). The largest available array is too small to follow the sign change and to decide between Raydec and array processing with certainty below 1.5 Hz. Between 2 and 4 Hz, all ellipticities are in agreement, including the corrected ellipticities calculated for the larger (PC5) and the smaller (PC2) arrays. Above 4 Hz, RTBF and RayDec are slightly diverging. The all-component techniques show a bump at 5 Hz which is not observed for the smaller array (Fig. S13b). Nonetheless, the values are compatible with the smaller array (PC2) up to 10 Hz. Above 10 Hz, the ellipticity tends to zero unlike the ellipticity of the smaller array PC2 located around −30° at 15 Hz.

In Fig. 17, the ellipticity correction and the incoherent noise indicators for ARDS are shown. The careful comparison between Figs 16(c) and 17(a) reveals that no significant modification is brought by the ellipticity correction, despite the curvature of $\delta$ between 0.8 and 1.8 Hz. At 1 Hz, $\delta$ is about 0.12 which corresponds to a correction of approximately 1° if $\xi_r$ is equal to 8°. For $|\delta|$ lower than 1, the ellipticity correction tends to 0 for a true ellipticity close to 0° or ±90° (Fig. 1). Between 3 and 4.5 Hz and between 6 and 10 Hz, the ellipticity correction is particularly efficient with a perfect fit of the corrected ellipticities obtained for arrays PC3 and PC2.

## 5.3 Cadarache, France

It is a quiet and flat site where limestones are outcropping (Garofalo et al. 2016a). Data sets labelled PC1 to PC4 (apertures from 30 to 270 m) are processed. Array PC3 which has an aperture of 270 m is analysed in details hereinafter. The results of the other data sets are available as supplementary materials in Fig. S15. There is no significant difference between the four methods whose results are shown in Fig. 18. A splitting into two curves is observed for all methods below 8 Hz for VBF and below 15 Hz for all 3C methods. Lateral variations which deviate from the 1-D assumption may explain this phenomenon, the details of which falls outside the scope of this work.

Raydec and RTBF ellipticities perfectly match between 2 and 50 Hz (Fig. 19), whilst ARTBF and ARDS ellipticities have a slightly lower absolute value below 10 Hz and a higher absolute value above 25 Hz. There is no significant difference between ARTBF and ARDS except at low frequency where ARTBF has higher absolute value than ARDS (−40° instead of −36° at 5 Hz). The two dispersion curves visible in Fig. 18 have slightly different ellipticity curves (about 4° higher in absolute value for the lowest velocity branch over the whole range, not shown). At high frequency (above 35 Hz in this case), the ellipticity tends to ±90° unlike the results in Figs 8, 13 and 16.
<figure><img src="image_17.png" /></figure>

**Figure 16.** Comparison of the ellipticity curves obtained simultaneously with the dispersion curves shown in Fig. 15: (a) RTBF, (b) ARTBF and (c) ARDS. The light blue curves are the mean (plain line) and one standard deviation (dashed lines) obtained for RayDec method over all stations of Grenoble array PC3 (Hobiger et al. 2009). A negative sign is added to the RayDec values for the sake of comparison. The purple curve found in all panels is manually picked from panel (c). The red and the dark blue curves are picked for the larger (PC5) and the smaller (PC2) arrays, respectively (Fig. S13). The colour scales are the same as in Fig. 8.

The ellipticity corrections shown in Fig. 20 moves the ellipticity towards lower absolute values except around 15 Hz where no modifications are observed. +5° and +12° are found at 6 and 30 Hz, respectively. At low frequency, the separation between ARDS and RayDec is increased by the correction. ARDS and ARTBF (not shown) ellipticities at low frequency (6 Hz) are found at approximately −33° which fits with the theoretical ellipticity of a half-space for a Poisson's ratio of 0.3, compatible with hard rocks. RayDec ellipticity of −43° at 6 Hz (or even lower for lower frequencies)

<figure><img src="image_18.png" /></figure>

**Figure 17.** ARDS results with corrected ellipticity for array PC3 (Grenoble): (a) ellipticity curves, (b) incoherent noise ratio for vertical components $\hat{R}_z/N$, (c) incoherent noise ratio for horizontal components $\hat{R}_h/N$ and (d) difference between $\hat{R}_z/N$ and $\hat{R}_h/N$ (δ). The light blue curves are the same as in Fig. 16. The purple curve is picked on the uncorrected results (Fig. 16) while the black curve is related to the corrected ellipticity. The red and the dark blue curves are picked for the larger (PC5) and the smaller (PC2) arrays, respectively (corrected ellipticities, Fig. S13). The colour scales are the same as in Figs 8 and 9.

corresponds to a Poisson's ratio lower than 0.1 which sounds unrealistic. A short discussion about the dependency of the ellipticity versus Poisson's ratio for a half-space can be found in Supporting Information (Section 1).

## 6 CONCLUSIONS

By questioning the structure of the incoherent noise assumed in Wathelet et al. (2018) and by considering the contributions of the incoherent noise to vertical and horizontal components separately, new formulations can be established for the ellipticity estimators $\hat{\xi}_h$, $\hat{\xi}_z$ and $\hat{\xi}_s$ ($\hat{e}_h$, $\hat{e}_z$ and $\hat{e}_s$ in the original work, linked to the new estimators by the tangent). In particular, the removal of the incoherent noise influence claimed by Wathelet et al. (2018) for $\hat{e}_s$ is disconfirmed. $\hat{\xi}_h$ and $\hat{\xi}_z$ corresponds to steering matrices where
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Comparison of the dispersion curves retrieved for Cadarache array PC3 (aperture = 270 m) with (a) VBF, (b) RTBF, (c) ARTBF and (d) ARDS. The black straight lines are the theoretical array limits, from left to right: $k_{min}/2$ (plain line), $k_{min}$ (dotted line), $\bar{k}_{max}/2$ (dotted line) and $\bar{k}_{max}$ (dashed line). The purple curve found in all panels is manually picked from panel (d). The blue curves are picked for the smaller array (PC4, Fig. S15).</figcaption></figure>

only the horizontal and the vertical components are modified in amplitude, respectively. The angular ellipticity used throughout this work suggests a third way where both components are modified in amplitude, leading to a new estimator, $\hat{\xi}_a$, which does not depend upon the incoherent levels like the others but upon the difference between the horizontal and the vertical noise ratios ($\delta$). None of the ellipticity estimators is noise-free. Hence, a correction strategy based on the computation of the beam power second derivative is proposed. The capabilities of a 3C cross-spectrum matrix (ARDS) are also compared to those of a projected two-component cross-spectrum matrix (ARTBF): only a minor improvement can be noted from a theoretical point of view.

The developed methods are tested on four different cases: a synthetic 1-D example and three experimental data sets collected during InterPACIFIC project (Mirandola, Grenoble and Cadarache). For the synthetic case, the ellipticity curves of the first six modes are well retrieved within the array wavenumber limits, in particular around 0° and ±90° where RTBF was not successful. Despite the weak theoretical background, ARDS (3N cross-spectral matrix) proves its superiority compared to ARTBF (2N cross-spectral matrix) whose results are partially contaminated by Love waves. The

ellipticity correction has not a significant effect for this noise-free data set. Adding incoherent noise on either the horizontal or the vertical components, the ellipticity results can be strongly biased even if a dispersion can still be delineated. In such a case, the ellipticity correction becomes an efficient tool to recover a correct ellipticity curve. For the experimental sites, realistic signed ellipticity curves are derived from the processing of several arrays of varying apertures. Love contamination of ARTBF results can be suspected only for Mirandola. For Grenoble, ARDS is less efficient than ARTBF in retrieving the first higher order mode. There are only minor differences between ARTBF and ARDS for Cadarache. The ellipticity correction always brings valuable and consistent modifications even if they are not always perfect.

The results are systematically compared to RayDec results computed on the same data sets. There is always an intermediate frequency range where RayDec results fit with the absolute ellipticity obtained with array processing. At low frequency, the single station and the array ellipticities always differ. For the synthetic case, the array ellipticity is more reliable. For the experimental sites, there is a body of evidence that the array ellipticity is preferable but it is impossible to prove it without additional studies.
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Comparison of the ellipticity curves obtained simultaneously with the dispersion curves shown in Fig. 18: (a) RTBF, (b) ARTBF and (c) ARDS. The light blue curves are the mean (plain line) and one standard deviation (dashed lines) obtained for RayDec method over all stations of Cadarache array *PC3* (Hobiger *et al.* 2009). A negative sign is added to the RayDec values for the sake of comparison. The purple curve found in all panels is manually picked from panel (c). The dark blue curves are picked for the smaller array (*PC4*, Fig. S15). The colour scales are the same as in Fig. 8.</figcaption></figure>

For the synthetic case, $\delta$ is flat and close to 0 at low frequency: there is no bias in the ellipticity determination. For all the experimental cases, this is never the case. $\delta$ is always deviating towards negative values even inside the array wavenumber limits. A negative $\delta$ corresponds to higher incoherent noise on the horizontal than on the vertical components. In a synthetic case, the sensors are perfectly oriented. This is not necessarily the case on the field where small compasses are used to orient sensors: an error of $\pm 5^\circ$ can be reasonably assumed. Further investigations should check the robustness of the proposed methods against orientation imperfections. At high frequency the ellipticity tends towards $0^\circ$ or $\pm 90^\circ$.

<figure><img src="image_21.png" /><figcaption>**Figure 20.** ARDS results with corrected ellipticity for array *PC3* (Cadarache): (a) ellipticity curves, (b) incoherent noise ratio for vertical components $R_z/N$, (c) incoherent noise ratio for horizontal components $R_h/N$ and (d) difference between $R_z/N$ and $R_h/N$ ($\delta$). The light blue curves are the same as in Fig. 19. The purple curves are picked on the uncorrected results (Fig. 19) while the black curve is related to the corrected ellipticity. The dark blue curve is picked for the smaller array (*PC4*, corrected ellipticities, Fig. S15). The colour scales are the same as in Figs 8 and 9.</figcaption></figure>

It is case and array dependent, in a currently unpredictable way. Understanding the role of the aliasing patterns can be one of the paths to follow in future studies.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

### suppl.data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## ACKNOWLEDGMENTS

This work would not have been possible without the numerous discussions I had with Bertrand Guillier during the preparation of this paper. I am also very grateful to Manuel Hobiger and Valerio
Poggy who reviewed the (long) paper and whose comments were
invaluable in improving the quality of this work. ISTerre is part of
Labex OSUG2020 (ANR10 LABX56).

DATA AVAILABILITY

The synthetic data set will be shared on reasonable request to the cor-
responding author. The experimental data underlying this article are
available at https://interpacific.geopsy.org. All tools developed in
this paper are distributed under the GNU Public License v3 and they
are available for download at https://www.geopsy.org (version>=
3.5.0).

REFERENCES

Asten, M.W. & Henstridge, J.D., 1984. Array estimators and the use of mi-
croseisms for reconnaissance of sedimentary basins, *Geophysics*, **49**(11),
1828–1837.

Bard, P.-Y., Bonnefoy-Claudet, S., Cornou, C., Fäh, D., Kristek, J. & Moczo, P., 2004. Simulation of seismic noise: report on parameter studies, *SESAME report D12.09*, Tech. rep., last downloaded from http://sesame.geopsy.org on March 2018.

Boaga, J., Cassiani, G., Strobbia, C.L. & Vignoli, G., 2013. Mode misiden-
tification in Rayleigh waves: ellipticity as a cause and a cure, *Geophysics*,
**78**(4), EN17–EN28.

Boore, D.M. & Toksöz, M.N., 1969. Rayleigh wave particle motion and crustal structure, *Bull. seism. Soc. Am.*, **59**(1), 331–346.

Capon, J., 1969. High-resolution frequency-wavenumber spectrum analysis, *Proc. IEEE*, **57**(8), 1408–1418.

Crameri, F., Shephard, G.E. & Heron, P.J., 2020. The misuse of colour in science communication, *Nat. Commun.*, **11**(1), 1–10.

Dunkin, J.W., 1965. Computation of modal solutions in layered, elastic media at high frequencies, *Bull. seism. Soc. Am.*, **55**(2), 335–358.

Fäh, D., Kind, F. & Giardini, D., 2001. A theoretical investigation of average H / V ratios, *Geophys. J. Int.*, **145**(2), 535–549.

Fäh, D., Stamm, G. & Havenith, H.-B., 2008. Analysis of three-component ambient vibration array measurements, *Geophys. J. Int.*, **172**(1), 199–213.

Fu, C.Y., 1946. Studies on seismic waves. II: Rayleigh waves in a superficial layer, *Geophysics*, **11**(1), 10–23.

Gal, M., Reading, A., Ellingsen, S., Koper, K., Burlacu, R. & Gibbons, S., 2016. Deconvolution enhanced direction of arrival estimation using one-and three-component seismic arrays applied to ocean induced microseisms, *Geophys. J. Int.*, **206**(1), 345–359.

Garofalo, F. et al., 2016a. InterPACIFIC project: comparison of invasive and non-invasive methods for seismic site characterization. Part I: intra-comparison of surface wave methods, *Soil Dyn. Earthq. Eng.*, **82**, 222–240.

Garofalo, F. et al., 2016b. InterPACIFIC project: comparison of inva-
sive and non-invasive methods for seismic site characterization. Part II:
inter-comparison between surface-wave and borehole methods, *Soil Dyn.
Earthq. Eng.*, **82**, 241–254.

Haskell, N.A., 1953. The dispersion of surface waves in multi-layered media, *Bull. seism. Soc. Am.*, **43**, 17–43.

Henderson, H.V. & Searle, S.R., 1981. On deriving the inverse of a sum of matrices, *SIAM Rev.*, **23**(1), 53–60.

Hobiger, M., Bard, P.-Y., Cornou, C. & Le Bihan, N., 2009. Single station determination of Rayleigh wave ellipticity by using the random decrement technique (RayDec), *Geophys. Res. Lett.*, **36**(L14303), 1–5.

Hobiger, M., Le Bihan, N., Cornou, C. & Bard, P.-Y., 2012. Multicomponent signal processing for Rayleigh wave ellipticity estimation: application to seismic hazard assessment, *IEEE Sig. Process. Mag.*, **29**(3), 29–39.

Hobiger, M. *et al.*, 2013. Ground structure imaging by inversions of Rayleigh wave ellipticity: sensitivity analysis and application to European strong-motion sites, *Geophys. J. Int.*, **192**(1), 207–229.

Hobiger, M. *et al.*, 2021. Site characterization of Swiss strong-motion sta-
tions: the benefit of advanced processing algorithms, *Bull. seism. Soc.
Am.*, **111**, 1713–1739.

Lachet, C. & Bard, P.-Y., 1994. Numerical and theoretical investigations on the possibilities and limitations of Nakamura's technique, *J. Phys. Earth*, **42**(5), 377-397.

Lacoss, R.T., Kelly, E.J. & Toksöz, M.N., 1969. Estimation of seismic noise structure using arrays, *Geophysics*, **34**(1), 21–38.

Lermo, J. & Chávez-García, F.J., 1994. Are microtremors useful in site response evaluation?, *Bull. seism. Soc. Am.*, **84**(5), 1350–1364.

Löer, K., Riahi, N. & Saenger, E.H., 2018. Three-component ambient noise beamforming in the Parkfield Area, *Geophys. J. Int.*, **213**(3), 1478–1491.

Love, A. E.H., 1911. *Some Problems of Geodynamics: Being an Essay to Which the Adams Prize in the University of Cambridge was Adjudged in 1911*, Univ. Press.

Malischewsky, P.G. & Scherbaum, F., 2004. Love's formula and H/V-ratio (ellipticity) of Rayleigh waves, *Wave Motion*, **40**, 57–67.

Maranò, S., Reller, C., Loeliger, H.-A. & Fäh, D., 2012. Seismic waves estimation and wavefield decomposition: application to ambient vibrations, *Geophys. J. Int.*, **191**(1), 175–188.

Maranò, S., Hobiger, M. & Fäh, D., 2017. Retrieval of Rayleigh wave ellipticity from ambient vibration recordings, *Geophys. J. Int.*, **209**(1), 334–352.

Molnar, S. et al., 2018. Application of microtremor horizontal-to-vertical spectral ratio (MHVSR) analysis for site characterization: state of the art, *Surv. Geophys.*, **39**, 613–631.

Mooney, H.M. & Bolt, B.A., 1966. Dispersive characteristics of the first three Rayleigh modes for a single surface layer, *Bull. seism. Soc. Am.*, **56**(1), 43–67.

Nakamura, Y., 1989. A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface, *Railw. Tech. Res. Inst., Quart. Rep.*, **30**(1), 25–33.

Nakamura, Y., 2000. Clear identification of fundamental idea of Nakamura's technique and its applications, in *Proceedings of the 12th World Conference On Earthquake Engineering*, Vol. **2656**, Auckland New Zealand.

Nocedal, J. & Wright, S., 2006. *Numerical Optimization*, Springer Science & Business Media.

Poggi, V. & Fäh, D., 2010. Estimating Rayleigh wave particle motion from three-component array analysis of ambient vibrations, *Geophys. J. Int.*, **180**(1), 251–267.

Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P., 2007. *Numerical Recipes: The Art of Scientific Computing*, Cambridge Univ. Press.

Riahi, N., Bokelmann, G., Sala, P. & Saenger, E.H., 2013. Time-lapse analysis of ambient surface wave anisotropy: a three-component array study above an underground gas storage, *J. geophys. Res.*, **118**(10), 5339–5351.

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**(1), 367–373.

Sánchez-Sesma, F.J. et al., 2011. A theory for microtremor H/V spectral ratio: application for a layered medium, *Geophys. J. Int.*, **186**(1), 221–225.

Scholte, J., 1947. On rayleigh waves in visco-elastic media, *Physica*, **13**(4), 245–250.

Sexton, J.L., Rudman, A. & Mead, J., 1977. Ellipticity of Rayleigh waves recorded in the Midwest, *Bull. seism. Soc. Am.*, **67**(2), 369–382.

Sezawa, K., 1927. Dispersion of elastic waves propagated on the surface of stratified bodies and on curved surfaces, *Bull. Earthq. Res. Inst., Tokyo Imper. Univ.*, **3**, 1–18.

Sezawa, K. & Kanai, K., 1938. Anomalous dispersion of elastic sur-
face waves. II, *Bull. Earthq. Res. Inst., Tokyo Imper. Univ.*, **16**,
683–689.

Stoneley, R., 1934. The transmission of Rayleigh waves in a heterogeneous medium, *Geophys. Suppl. Mon. Not. Astron. Soc.*, **3**(6), 222-232.

Strutt, J.W., 1885. On waves propagated along the plane surface of an elastic solid, *Proc. Lond. Math. Soc.*, **s1-17**(1), 4–11.

Takagi, R., Nishida, K., Maeda, T. & Obara, K., 2018. Ambient seismic noise wavefield in Japan characterized by polarization analysis of Hi-net records, *Geophys. J. Int.*, **215**(3), 1682–1699.

Tanimoto, T., Ishimaru, S. & Alvizuri, C., 2006. Seasonality in particle motion of microseisms, *Geophys. J. Int.*, **166**(1), 253–266.

Thomson, W.T., 1950. Transmission of elastic waves through a stratified solid medium, *J. appl. Phys.*, **21**(2), 89-93.
Wagner, G.S., 1996. Resolving diversely polarized, superimposed signals in three-component seismic array data, *Geophys. Res. Lett.*, **23**(14), 1837–1840.

Wathelet, M., Jongmans, D., Ohrnberger, M. & Bonnefoy-Claudet, S., 2008. Array performances for ambient vibrations on a shallow structure and consequences over V s inversion, *J. Seismol.*, **12**(1), 1–19.

Wathelet, M., Guillier, B., Cornou, C., Roux, P. & Ohrnberger, M., 2018. Rayleigh wave three-component beamforming: signed ellipticity assessment from high-resolution frequency-wavenumber processing of ambient vibration arrays, *Geophys. J. Int.*, **215**(1), 507–523.

## APPENDIX A: DERIVATIVES OF BEAM POWER VERSUS ELLIPTICITY

If we define $E'_{2,h} = (-j\mathbf{q} \mathbf{0})^T$ and $E'_{2,z} = (\mathbf{0} \mathbf{q})^T$, we have from eq. (47) and Table 1,

$$ 
\begin{align} 
\frac{\partial P_{2,h}(\vec{k}, \xi)}{\partial \xi} &= -2 P_{2,h}^2 \frac{1}{\cos^2 \xi} \Re(E'_{2,h} F_p^{-1} E_{2,h}), \nonumber \\ 
\frac{\partial^2 P_{2,h}(\vec{k}, \xi)}{\partial \xi^2} &= \frac{2}{P_{2,h}} \left(\frac{\partial P_{2,h}}{\partial \xi}\right)^2 + \tan \xi \frac{\partial P_{2,h}}{\partial \xi} - 2 P_{2,h}^2 \frac{1}{\cos^4 \xi} (E'_{2,h} F_p^{-1} E'_{2,h})^\dagger, \tag{A1} \\ 
\frac{\partial P_{2,z}(\vec{k}, \xi)}{\partial \xi} &= 2 P_{2,z}^2 \frac{1}{\sin^2 \xi} \Re(E'_{2,z} F_p^{-1} E_{2,z}), \nonumber \\ 
\frac{\partial^2 P_{2,z}(\vec{k}, \xi)}{\partial \xi^2} &= \frac{2}{P_{2,z}} \left(\frac{\partial P_{2,z}}{\partial \xi}\right)^2 - 2 \cot \xi \frac{\partial P_{2,z}}{\partial \xi} - 2 P_{2,z}^2 \frac{1}{\sin^4 \xi} (E'_{2,z} F_p^{-1} E'_{2,z})^\dagger. \nonumber 
\end{align} 
$$

These formulas are also valid for $P_{3,h}$ and $P_{3,z}$ (by replacing $\mathbf{q}$ by $\mathbf{V}$ in the steering matrix). The derivatives of $P_{2,s}$ can be formulated from the values deduced for $P_{2,h}$.

$$ 
\begin{align} 
\frac{\partial P_{2,s}(\vec{k}, \xi)}{\partial \xi} &= 2 \frac{\tan \xi}{\cos^2 \xi} P_{2,h}^2 + 2 \tan^2 \xi P_{2,h} \frac{\partial P_{2,h}}{\partial \xi} \dagger \nonumber \\ 
\frac{\partial^2 P_{2,s}(\vec{k}, \xi)}{\partial \xi^2} &= 2 \left[ \frac{1}{\cos^2 \xi} \left( \frac{1}{\cos^2 \xi} + 2 \tan^2 \xi \right) P_{2,h}^2 + 4 \frac{\tan \xi}{\cos^2 \xi} P_{2,h} \frac{\partial P_{2,h}}{\partial \xi} + \tan^2 \xi \left( \left[ \frac{\partial P_{2,h}}{\partial \xi} \right]^2 + P_{2,h} \frac{\partial^2 P_{2,h}}{\partial \xi^2} \right) \right] \dagger \tag{A2} 
\end{align} 
$$

The analytical expression of the derivatives of $P_{2,h}$, $P_{2,z}$ and $P_{2,s}$ are slightly more complex than for $P_{3,a}$ or $P_{2,a}$. However, the cost is not significantly larger because there are only a few additional real operations, with no proportionality to $N$.

The iterative search can be started from 0° for $P_{2,a}$, $P_{3,a}$, $P_{2,h}$ and $P_{3,h}$. For the other beam powers, the first derivative is null or infinite at $\xi = 0$. A starting value of 45° is chosen. As shown in Fig. 2(c), the curve of $P_{2,s}(\xi)$ has always two maxima located at exactly the same absolute values but with opposite signs (eq. 35). Of the two maxima, the one with the highest amplitude corresponds to the true ellipticity. For the next optimizations, the previous optimum value is reused. In a practically unpredictable way, the obtained ellipticity can be positive or negative according to the path followed during the search. The sign of the obtained optimum ellipticity is inverted and a second optimization is run which must provide the same absolute value as the first result. Generally, the iterative process returns immediately because the required precision is directly achieved. From the two optimum values, only the ellipticity for which $P_{2,s}$ has the highest amplitude is kept.

## APPENDIX B: DERIVATIVES OF BEAM POWERS VERSUS THE RADIAL WAVENUMBER

The beamforming array response function is by definition

$$ B(\vec{k}) = \frac{1}{N} \sum_{n=1}^{N} q_n(\vec{k}). \tag{B1} $$

Hence the beam pattern is

$$ |B(\vec{k})|^2 = B^* B = \frac{1}{N^2} \sum_{n=1}^{N} \sum_{m=1}^{N} q_n^*(\vec{k}) q_m(\vec{k}). \tag{B2} $$

This is the mean of the element of the matrix $qq^*$ which can be written with matrix notation

$$ |B(\vec{k})|^2 = \frac{1}{N^2} \mathbf{1}^* \mathbf{q} \mathbf{q}^* \mathbf{1}, \tag{B3} $$

where $\mathbf{1}$ is a $N \times 1$ matrix filled with 1. If $k = |\vec{k}|$, the beamforming array response pattern derivatives along the wavenumber radial coordinate in a polar coordinate system are

$$ 
\begin{align} 
\frac{\partial}{\partial k} |B(\vec{k})|^2 &= \frac{2}{N^2} \Re(\mathbf{1}^* \mathbf{q} \frac{\partial \mathbf{q}^*}{\partial k} \mathbf{1}), \nonumber \\ 
\frac{\partial^2}{\partial k^2} |B(\vec{k})|^2 &= \frac{2}{N^2} \Re(\mathbf{1}^* \mathbf{q} \frac{\partial^2 \mathbf{q}^*}{\partial k^2} \mathbf{1} + \mathbf{1}^* \frac{\partial \mathbf{q}}{\partial k} \frac{\partial \mathbf{q}^*}{\partial k}). \tag{B4} 
\end{align} 
$$
The derivative of $q^*$ versus $k$ is

$$ \frac{\partial q^*}{\partial k} = j r_1^* \odot q^* \dagger, \qquad (B5) $$

where $r_1$ is a $N \times 1$ matrix containing the projections of the sensor coordinates on the direction $\theta$. The operator $\odot$ denotes the Hadamard product of two matrices (or the element-wise product). The first two derivatives of the array response pattern become

$$ \begin{aligned} \frac{\partial}{\partial k} |B(\vec{k})|^2 &= \frac{-2}{N^2} \mathbf{1}^* \mathfrak{S}(\mathbf{q}(r_1^* \odot \mathbf{q}^*))\mathbf{1}, \\ \frac{\partial^2}{\partial k^2} |B(\vec{k})|^2 &= \frac{-2}{N^2} \mathbf{1}^* \mathfrak{S}(\mathbf{q}(r_1^* \odot r_1^* \odot \mathbf{q}^*) - (r_1 \odot \mathbf{q})(r_1^* \odot \mathbf{q}^*))\mathbf{1}. \end{aligned} \qquad (B6) $$

At the origin ($\vec{k} = \vec{0}$), the first derivative is null and the second derivative reduces to

$$ \frac{\partial^2}{\partial k^2} |B(\vec{k})|^2(\vec{0}) = \frac{2}{N^2} [\mathbf{1}^* \mathbf{r}_1 \mathbf{r}_1^* \mathbf{1} - N(\mathbf{r}_1^* \odot \mathbf{r}_1^*) \mathbf{1}]. \qquad (B7) $$

The beam power derivatives along the wavenumber radial coordinate in the same polar coordinate system are similar to the derivative versus $\xi$ developed in eq. (47), that is,

$$ \begin{aligned} \frac{\partial P}{\partial k} &= -2P^2 \mathfrak{S} \left( \frac{\partial E^*}{\partial k} F^{-1} E \right) \quad \text{and} \\ \frac{\partial^2 P}{\partial k^2} &= \frac{2}{P} \left[ \frac{\partial P}{\partial k} \right]^2 - 2P^2 \mathfrak{S} \left( \frac{\partial^2 E^*}{\partial k^2} F^{-1} E + \frac{\partial E^*}{\partial k} F^{-1} \frac{\partial E}{\partial k} \right)^{\dagger}, \end{aligned} \qquad (B8) $$

where

$$ \frac{\partial E^*}{\partial k} = j r^* \odot E^*; \qquad (B9) $$

$r$ is either $r_2$ ($2N \times 1$) or $r_3$ ($3N \times 1$), containing two or three times the same $N$ projections of the sensor coordinates on direction $\theta$. $P$, $F$ and $E$ are generic quantities. Eq. (B8) applies to any of the beam powers proposed in this work. The first two derivatives of the beam power become

$$ \begin{aligned} \frac{\partial P}{\partial k} &= 2P^2 \mathfrak{S}((r^* \odot E^*)F^{-1}E), \\ \frac{\partial^2 P}{\partial k^2} &= \frac{2}{P} \left[ \frac{\partial P}{\partial k} \right]^2 + 2P^2 \mathfrak{S}((r^* \odot r^* \odot E^*)F^{-1}E) - 2P^2(r^* \odot E^*)F^{-1}(r \odot E). \end{aligned} \qquad (B10) $$

Like the derivative of the beam power versus the ellipticity, the evaluation of the beam power and its two derivatives costs a bit more than twice the computation of the beam power value alone.