\makeatletter
\@ifundefined{HCode}
{\documentclass[screen, CRGEOS, Unicode, Thematic]{cedram}
\usepackage[square,comma,authoryear]{natbib}
\renewcommand\bibsection{\section*{\refname}}
\newenvironment{noXML}{}{}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{}%\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.7\linewidth}{#1}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\ndash{\text{--}}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\def\0{\phantom{0}}
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
\RequirePackage{etoolbox}
\usepackage{upgreek}
\def\jobid{crgeos20221200}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\skip\footins=16pt
\newcounter{runlevel}
\usepackage{multirow}
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\let\ubreak\break
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\def\back#1{}
\def\rmmu{\upmu}
\DOI{10.5802/crgeos.254}
\datereceived{2022-12-08}
\daterevised{2023-12-05}
\datererevised{2024-02-09}
\dateaccepted{2024-01-29}
\ItHasTeXPublished
}
{
\PassOptionsToPackage{authoryear}{natbib}
\documentclass[crgeos]{article}
\usepackage[T1]{fontenc}
\def\CDRdoi{10.5802/crgeos.254}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\let\ubreak\relax
\let\newline\relax
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
}
\makeatother

\usepackage{upgreek}

\dateposted{2024-02-22}
\begin{document}

\begin{noXML}

%\makeatletter
%\def\TITREspecial{\relax}
%\def\cdr@specialtitle@english{New Developments in Passive Seismic Imaging and Monitoring}
%%\def\cdr@specialtitle@french{Nouveaux d\'eveloppements en mati\`ere d'imagerie sismique passive et de surveillance}
%\makeatother

\title{Investigating the lateral resolution of the Rayleigh wave focal
spot imaging technique using two-dimensional acoustic simulations}

\alttitle{\'Etude de la r\'esolution lat\'erale de l'imagerie par tache
focale des ondes de Rayleigh \`a l'aide de simulations acoustiques
bidimensionnelles}

\author{\firstname{Bruno} \lastname{Giammarinaro}\CDRorcid{0000-0001-9968-1915}\IsCorresp}
\address{Institute of Seismology, University of Helsinki, Helsinki,
Finland}
\curraddr[B. Giammarinaro]{LabTAU, INSERM, Centre L\'eon B\'erard,
Universit\'e Lyon~1, Lyon, France}
\email[B. Giammarinaro]{bruno.giammarinaro@inserm.fr}

\author{\firstname{Christina} \lastname{Tsarsitalidou}\CDRorcid{0000-0003-1404-4648}}
\addressSameAs{1}{Institute of Seismology, University of Helsinki,
Helsinki, Finland}
\email[C. Tsarsitalidou]{christina.tsarsitalidou@helsinki.fi}

\author{\firstname{Gregor} \lastname{Hillers}\CDRorcid{0000-0003-2341-1892}}
\addressSameAs{1}{Institute of Seismology, University of Helsinki,
Helsinki, Finland}
\email[G. Hillers]{gregor.hillers@helsinki.fi}

\shortrunauthors

\keywords{\kwd{Numerical modeling}\kwd{Wave propagation}\kwd{Seismic
noise}\kwd{Cavity}\kwd{Focal spot}\kwd{Resolution}\kwd{Seismic
imaging}}

\altkeywords{\kwd{Mod\'elisation num\'erique}\kwd{Propagation des
ondes}\kwd{Bruit sismique}\kwd{Cavit\'e}\kwd{Tache
focal}\kwd{R\'esolution}\kwd{Imagerie sismique}}

\thanks{Academy of Finland (decision number 322421)}

\begin{abstract}
We investigate the lateral resolution power of the seismic Rayleigh
wave focal spot imaging technique. We use two-dimensional acoustics
simulations in a closed cavity for the passive Green's function and
focal spot reconstruction. Four different velocity distributions target
different resolution aspects. The finite data range that is necessary
to constrain the Bessel function model controls the lateral spreading
of material contrasts, the distinction of two objects on sub-wavelength
scales, and the image quality of complex random media. Good data
quality from dense networks supports short range estimates and
super-resolution.
\end{abstract}

\begin{altabstract}
Nous \'etudions la r\'esolution lat\'erale de la technique d'imagerie
utilisant les taches focales des ondes sismiques de Rayleigh. Nous
utilisons des simulations acoustiques bidimensionnelles dans une
cavit\'e ferm\'ee pour la reconstruction passive taches focales et des
fonctions de Green. Quatre distributions de vitesse diff\'erentes
ciblent diff\'erents aspects de la r\'esolution. La plage finie de
donn\'ees n\'ecessaire pour contraindre le mod\`ele de fonction de
Bessel contr\^ole l'\'etalement lat\'eral des contrastes observ\'es
dans mat\'eriels, la distinction de deux objets sur des \'echelles
inf\'erieures \`a la longueur d'onde, et la qualit\'e d'image dans des
milieux al\'eatoires complexes. La bonne qualit\'e des donn\'ees issue
des r\'eseaux denses autorise les estimations \`a courte port\'ee et la
super-r\'esolution.
\end{altabstract}

\maketitle
\vspace*{.35pc}

\twocolumngrid

\end{noXML}

\section{Introduction}

Disadvantages of earthquake tomography associated with limited
illumination can now be compensated by ambient noise tomography with
its flexible virtual source and receiver configurations. Both
approaches invert far-field observations of travel time differences,
obtained from earthquake seismograms \citep{Tromp2005, LiuQ2012} or
from passive Green's function reconstructions \citep{Shapiro2004,
Sabra2005}, for a model of the velocity structure.

Modern dense seismic arrays support alternative local surface wave
speed estimations from noise correlation functions \citep{Lin2009,
LinFC2011}, which includes the large scale application of the frequency
domain spatial autocorrelation (SPAC) method \citep{Aki1957,
Ekstrom2009, Ekstrom2014} that is otherwise typically applied to local
sparse array data \citep{Asten2006}.
Dense arrays can now contain on the order of 1000 sensors, which
facilitates the proper sampling of the noise correlation amplitude
distribution at near-field distances. At zero lag time, the time domain
representation of the spatial autocorrelation field is referred to as
focal spot, which contains the same information as SPAC and can be
analyzed using the same mathematical tools \citep{Cox1973, Yokoi2008,
Tsai2010, Haney2012, Haney2014}.

Focal spot analysis is used in applications that work with a high
sensor density. Focal spots have first been studied in time-reversal
experiments in underwater acoustics and medical imaging
\citep{Fink1997}. Noise correlation medical imaging transferred the
noise seismology approaches to passive elastography, first using
far-field surface waves along muscle fibers \citep{Sabra2007}, later
using refocusing shear waves. Different properties of the shear wave
focal spot have been analyzed including its cross-section width
\citep{Catheline2008, Gallot2011}, two-dimensional shape
\citep{Benech2013, Brum2015}, and its tip curvature
\citep{Catheline2013}, which can all be reconstructed using MRI
\citep{Zorgani2015} and ultrasound \citep{Barrere2020} speckle tracking
methods. Importantly, \citet{Zemzemi2020} demonstrated that the ability
to discriminate two objects is not controlled and hence limited by the
shear wavelength, but instead by the ultrasonic frequency and the pixel
density. In seismology, this corresponds to the array station density.

\citet{Hillers2016} first applied Rayleigh wave focal spot imaging in
seismology to image lateral velocity variations in a fault zone
environment. As for SPAC, the phase velocity is estimated from the
focal spot shape using Bessel function models. Using numerical
time-reversal experiments based on a Green's function calculator for
one-dimensional layered media \citep{Cotton1997},
\citet{Giammarinaro2023} demonstrated the feasibility to accurately
estimate phase velocity and dispersion from noisy focal spots. Results
for non-isotropic surface wave illumination showed that the wave speed
bias is negligible if the sensors are isotropically distributed,
compatible with the SPAC results by \citet{Nakahara2006}. Moreover, the
effect of interfering $P$~waves can be mitigated using a data or
fitting range around one wavelength. Whereas these results from
\citet{Giammarinaro2023} suggest the overall robustness and utility of
the focal spot method for seismic imaging applications, most notably
because of the increase in depth resolution, the study could not
address lateral resolution. Considering that super-resolution can be
obtained with tip-curvature measurements \citep{Zemzemi2020}, it is
important to assess the effect of the data range on speed estimates
from densely sampled seismic focal spots.

Here we study systematically the lateral focal spot resolution using
numerical experiments. We perform two-dimensional acoustics simulations
to reconstruct the Green's function from reverberating wave fields
(Section~\ref{sec:method}). The ambient field generated in a chaotic
closed cavity \citep{Draeger1999} yields results that are equivalent to
results from open media noise correlation \citep{Derode2003}. The
obtained Green's function is identical and is here therefore taken as a
proxy for seismic vertical--vertical component Rayleigh wave
correlations \citep{Sanchez-Sesma2006, Haney2012}. We work with a
constant number of grid points and a fixed reference frequency. We
implement four test cases and vary the data range or fitting distance
$r_{\mathrm{fit}}$ to investigate the effect on the resolution of the
velocity structure. These cases include a homogeneous control
experiment (Section~\ref{sec:res_homo}), an interface between two
half-spaces (Section~\ref{sec:res_halfspaces}), circular inclusions
(Section~\ref{sec:res_inclusions}), and heterogeneous or random
velocity distributions (Section~\ref{sec:res_random}). In
Section~\ref{sec:discussion} we discuss the different resolution
aspects that are investigated with the variable configurations for a
comprehensive evaluation of the seismic Rayleigh wave focal spot
imaging performance.

\section{Method} \label{sec:method}
\subsection{Synthetic experiments}

This study is based on synthetic diffuse wave fields generated in a
closed cavity. The resulting correlation functions are equivalent to
noise correlations in open media \citep{Derode2003}. Simulations are
performed using the function \textit{kspaceFirstOrder2D} from the
MATLAB toolbox \textit{kWave}\ \citep{Treeby2018}. This function solves
a system of first-order acoustic equations for the conservation of mass
and momentum using a wavenumber $k$-space pseudospectral method. The
two-dimensional medium is composed of $500 \times 500$ grid points that
are spaced in the $x$ and $y$ direction by $\mathrm{d}x = \mathrm{d}y =
0.1$~km (Figure~\ref{fig:cavity}a). The background wave speed in the
cavity is $V_0 = 2$~km/s. The closed cavity is implemented using
different densities outside ($\rho_{\mathrm{out}} = 59~\text{kg/m}^3$)
and inside ($\rho_{\mathrm{in}} = 2950~\text{kg/m}^3$) the cavity.
Choosing $\rho_{\mathrm{out}}$ to be 2\% of $\rho_{\mathrm{in}}$
creates strongly reflecting \mbox{boundaries} from impedance contrast.
It supports a homogeneous stability regime of the \textit{kWave}
simulations through a constant Courant--Friedrichs--Lewy number. This
number depends on the wave speed, hence varying only the density
mitigates potentially problematic stability conditions associated with
a large change in wave speed. We cannot exclude the occurrence of weak
numerical dispersion in some of the heterogeneous case studies, but the
overall consistency of the synthesized Green's functions and focal
spots in the different experiments suggests that this effect does not
govern our results. Inside the cavity we select a square target domain
consisting of $151 \times 151$ grid points where we record the
solution. In this region we define the different velocity distributions
introduced in Section~\ref{sec:AcMedium}. Results for each of the four
cases discussed in Section~\ref{sec:Results} are obtained by averaging
11 independent wave field simulations. Each simulation starts with a
point source at a different position inside the cavity 
(Figure~\ref{fig:cavity}a). The source term is defined as a
time-varying mass source (Figures~\ref{fig:cavity}b, c)
emitting a 1~s long pulse. The pressure is a differentiation of this
signal with the frequency range centered at 1~Hz. The wave field is
recorded for 300~s inside the square sub-domain with a 100~Hz sampling
frequency.\vspace*{-.4pc}

\begin{figure*}
\includegraphics{fig01}
{\vspace*{-.15pc}}
\caption{\label{fig:cavity}Configuration of the numerical experiments.
(a)~Representation of the closed cavity with a background velocity $V_0
= 2$~km/s. The gray area represents the outer part of the medium with
the same velocity but with an impedance contrast to trap the waves in
the cavity. The red line is the boundary of the closed cavity. The
black square indicates the target area where the results are recorded.
Inside the cavity the color corresponds to the input wave speed. Red
dots indicate the source positions for the 11 realizations. (b)~Times
series and (c)~normalized power spectrum of the emitted pulse used for
each simulation. (d)~Example of a full time-series recorded in the
\mbox{cavity}.}
{\vspace*{-.42pc}}
\end{figure*}

\subsection{The acoustic medium in the target domain}
\label{sec:AcMedium}

The acoustic wave speed distribution $V(\mathbf{x})$ is defined as a
spatial function in the square target domain
{\begin{equation}
V(\mathbf{x})= V_0 (1 + \xi(\mathbf{x})), \label{eq:Cmap}
\end{equation}}\unskip
where $V_0 = 2$~km/s is the background wave speed in the cavity,
$\mathbf{x}$ is the position, and $\xi$ is the relative change in wave
speed, i.e., the parameter that controls the medium heterogeneity. In
the following sections, the 2~km background wavelength at 1~Hz is
denoted~$\lambda_0$. In Section~\ref{sec:Results} we begin with a
control experiment of a homogeneous medium with $\xi = 0$ to study the
overall system response.\vspace*{-.4pc}

\subsection{Lateral spreading across two welded{\hfill\break}
half-spaces}

We modify the homogeneous control experiment and replace the 2~km/s
velocity in the left half of the target domain by an increased 2.2~km/s
value. This creates two half-spaces and allows us to investigate the
lateral resolution as the imaging method induces spreading or averaging
across the sharp \mbox{interface}. Such sharp lateral velocity
contrasts can occur across bimaterial interfaces in fault zone
environments \citep{weertman_unstable_1980, ben-zion_response_1989}, or
in the contact region of intrusions with the host rock
\citep{Chamarczuk2019}.
\vspace*{-.25pc}

\subsection{Resolution of circular inclusions}

Next we perform a classic resolution test and study the power of the
method to separate two individual entities in an image. For this we
impose three pairs of circular inclusions separated by $0.25\lambda_0$,
$0.5\lambda_0$ and $1\lambda_0$. The inclusions have a diameter of
$1\lambda_0$. We test two different sets, where one set of inclusions
is stiffer than the background with $\xi = 25\%$, and the other set has
more compliant inclusions compared to the background with $\xi =
-25\%$. Such a ``two-body problem'' is typically studied in gel
phantom experiments performed in medical ultrasound imaging or medical
imaging for tumor detection \citep{Catheline2013, Zemzemi2020}. It is a
less common configuration in seismology where so-called checkerboard
tests are typically employed to quantify the resolution of a tomography
configuration. Circular cross-section features can occur in the context
of magmatic intrusions or conduits. However, as said, with this
configuration we can study the lateral resolution defined as the
minimum distance between objects that the method allows to discriminate
in an image. This and the other cases are studied using different data
ranges $r_{\mathrm{fit}}$ discussed in Section~\ref{sec:Processing}.
\vspace*{-.25pc}

\subsection{Randomly distributed wave speeds}

Last we consider random media using a functional form that is often
used to parameterize the heterogeneous distributions of variables such
as wave speed, stress, or frictional properties in Earth materials
\citep{Frankel1986, Holliger1992, Mai2002, Ripperger2007, Hillers2007,
Sato2012, Obermann2016}. We define $\xi(\mathbf{x})$ through a spatial
2D inverse Fourier transform as
{\begin{equation}
\xi(\mathbf{x})=\text{FT}^{-1}\left[\sqrt{P(\mathbf{k})}
\mathrm{e}^{\mathrm{i}\phi(\mathbf{k})}\right], \label{eq:Xi}
\end{equation}}\unskip
where $\mathbf{k}$ is the spatial wavenumber of the 2D distribution,
$P(\mathbf{k})$ is the power spectral density, and $\phi(\mathbf{k})$
is a random distribution of the phase between~0 and~$2\uppi$. The
random distribution is calculated using the Python \textit{randn}
function coupled with a seed fixed to~3. Fixing the seed allows
to randomly generate the phase and to keep the same distribution to
observe the effect of the control parameters. The power spectral
density $P(\mathbf{k})$ follows the von~Karman probability function
\citep{Sato2012}
{\begin{equation}
P(\mathbf{k}) = \frac{4\uppi\Gamma[\kappa+1] \varepsilon^2 a^2}
{\Gamma[\kappa](1+a^2 \|\mathbf{k}\|^2)^{\kappa+1}}. \label{eq:P}
\end{equation}}\unskip
The correlation length of the modeled parameter is $a$, $\varepsilon$
governs the contrast in the medium, $\kappa$ defines the sharpness of
the spectral decay, and $\Gamma$ is the Gamma function. We use
Equation~(\ref{eq:P}) to generate nine different media with variable
$a$ and $\kappa$ (Table~\ref{tab:RandomParameters}), and constant
$\varepsilon = 150~\text{km}^{-1}$.
{\vspace*{-.35pc}}

\begin{table*}[t!]%tab1
\tabcolsep=8.8pt
\caption{\label{tab:RandomParameters}Parameters for the generation of
different heterogeneous media using Equation~(\ref{eq:P})}
\begin{tabular}{cccccccccc}
\thead
Medium & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9\\
\endthead
$a$ (km) & 1 & 1 & 1 & 2 & 2 & 2 & 3 & 3 & 3\\
$\kappa$ & 0.1 & 0.3 & 0.6 & 0.1 & 0.3 & 0.6 & 0.1 & 0.3 & 0.6
\botline
\end{tabular}
\tabnote{Where $a$ is the correlation length and $\kappa$ governs the
sharpness of the spectral decay. The corresponding wave speed
distributions are displayed in the top row of
Figure~\ref{fig:V_estimation}.}
\end{table*}

\subsection{Data processing and wave speed estimation}
\label{sec:Processing} \label{sec:data_proc}

Data processing is performed in a Python3.8 environment. For each
simulation, we analyze the diffusive parts of the wave field
(Figure~\ref{fig:Example_proc}a) by focusing on simulated time series
data between 100~s and 300~s (Figure~\ref{fig:cavity}d). Considering
the distances \mbox{between} the sources and chaotic cavity borders,
this corresponds to a range between 10 and 30 mean-free times. We
follow \citet{Giammarinaro2023} and filter the traces with a Gaussian
filter centered at $1$~Hz and with a width of 3.2\% to estimate the
central frequency dependent phase velocity. In the limit, the analysis
could be performed for a harmonic case, however, the narrow filter
stabilizes the results. As demonstrated by \citet{Giammarinaro2023}, a
wider frequency-band filter leads to group velocity estimates, and it
biases the focal spot reconstruction by averaging over frequencies
which leads to an apparent attenuation. We compute the normalized
cross-correlation between each sensor pair to extract the spatial
autocorrelation amplitude fields at zero lag time. Results are stacked
over the 11 realizations for each case. This yields at short distances
around a reference station the large-amplitude feature referred to as
focal spot (Figure~\ref{fig:Example_proc}c). The shape is linked to the
wave speed through the imaginary part of the Green's function. From
each focal spot we estimate the free parameter wave speed $V =
\omega/k$ by fitting the azimuthal averaged and normalized data to the
$J_0(kr)$ model using a nonlinear least squares regression algorithm
\citep{Hillers2016, Giammarinaro2022Data, Giammarinaro2023}. Again, the
2D acoustic configuration yields results that are equivalent to the
lateral propagation of Rayleigh surface waves. The $J_0(kr)$ model
equally describes the vertical--vertical component of the Rayleigh wave
focal spot \citep{Haney2012, Haney2014}, and $V$ is thus equivalent to
the Rayleigh wave phase velocity~$c_R$.

\begin{figure*}
\includegraphics{fig02}
\caption{\label{fig:Example_proc}Workflow stage examples for the
homogeneous medium in the target area. (a)~Snapshot of a diffuse wave
field. (b)~Time-reversed space-time wave field obtained by
cross-correlation. (c)~Spatial autocorrelation with the focal spot at
the center obtained at 1~Hz after stacking over 11 realizations.
(d)~Simulated 1~Hz focal spot data (gray) and nonlinear regression
results for $r_{\mathrm{fit}} = 0.5\lambda_0$ (orange) and
$r_{\mathrm{fit}} = 1\lambda_0$ (red). In this and all subsequent
figures the unit $\lambda$ is equal to the 1~Hz wavelength $\lambda_0$
for the reference $V_0 = 2$~km/s.}
{\vspace*{-.5pc}}
\end{figure*}

This process is performed for the 22,801 focal spots, and each obtained
$V$ estimate is associated with the location of the reference station.
This imaging concept thus compiles velocity distributions across dense
arrays without solving a tomography inverse problem. Importantly, we
choose two different data ranges $r_{\mathrm{fit}}$ of 1~km and 2~km
associated with $0.5\lambda_0$ and $1\lambda_0$ at 1~Hz
(Figure~\ref{fig:Example_proc}d). Away from edges, this corresponds to
80 and 314 samples, respectively. We vary $r_{\mathrm{fit}}$ because it
is a critical tuning parameter, and values around one wavelength yield
overall stable results \citep{Giammarinaro2023}. Larger values can
stabilize a regression for noisy signals, but the short distance focal
spot imaging concept essentially invites the limitation of
$r_{\mathrm{fit}}$ for improved resolution. This refers to improved
depth resolution as $c_R$ values can be estimated at wavelengths that
cannot be studied with tomography \citep{TsarsitalidouEGU2021,
Giammarinaro2023}, but also to the lateral resolution investigated
here.

\subsection{Error estimation}

An important advantage of the focal spot approach is the local
assessment of the measurement uncertainty $\epsilon$ that is estimated
following \citet{Aster2019book}
{\begin{equation}
\epsilon_k = \sqrt{\frac{\mathit{RSS}}{N-2} C_k}, \label{eq:error_k}
\end{equation}}\unskip
where \textit{RSS} is the residual sum of squares, $N$ is the number of
samples, and $C_k$ is the diagonal element of the parameter covariance
matrix $C$ associated with the wavenumber~$k$. Rules of error
propagation apply to yield the wave speed error $\epsilon_V$
{\begin{equation}
\epsilon_V = V\epsilon_k / k. \label{eq:error_V}
\end{equation}}\unskip
In addition to this estimate we can compare imaged with ground truth
input values to discuss resolution.

\section{Results} \label{sec:Results}
\subsection{The homogeneous reference case} \label{sec:res_homo}

The first test is the homogeneous control experiment which illustrates
basic features of the approach. The simulations yield diffuse wave
fields that can be used for focal spot imaging. 
Figure~\ref{fig:Example_proc}(a) shows a snapshot of a diffuse wave
field, Figure~\ref{fig:Example_proc}(b) a time-space representation of
the time-reversed correlation wave field,
Figure~\ref{fig:Example_proc}(c) shows a focal spot of one realization,
and Figure~\ref{fig:Example_proc}(d) shows the results of the nonlinear
regression following the data processing described in
Section~\ref{sec:data_proc}. In Figure~\ref{fig:Example_proc}(b) we can
see the refocusing and diverging waves of the Green's function after
cross-correlation. However, the reconstruction is not perfect which is
indicated by the small-amplitude fluctuations. The time domain
autocorrelation field shows the focal spot at small distances around
the origin (Figure~\ref{fig:Example_proc}c). The irregularity of the
white zero-crossing contours again illustrates fluctuations and
imperfect reconstruction after stacking over 11 realizations. We
attribute this to non-perfectly diffuse wave fields associated with the
modes in the cavity. However, studies of non-diffuse wave fields
\citep{Nakahara2006, Giammarinaro2023} \mbox{demonstrate} that
isotropically distributed sensors yield unbiased wave speed estimates.
This condition is not met along the boundaries of the domain which
yields larger fluctuations in the estimates.
Figure~\ref{fig:Example_proc}(d) displays results from the nonlinear
regression using the two data distances $r_{\mathrm{fit}} =
0.5\lambda_0$ and $r_{\mathrm{fit}} = 1\lambda_0$. The spread in the
data with increasing distance corresponds to the fluctuations from
Figure~\ref{fig:Example_proc}(c).

The 2D velocity distributions obtained with $r_{\mathrm{fit}} =
0.5\lambda_0$ and $r_{\mathrm{fit}} = 1\lambda_0$ are displayed in
Figures~\ref{fig:C_control}(a) and~(b),
respectively. The key feature in these images are the speckle patterns.
They illustrate that the imperfect reconstruction affects the
measurement process to produce these fluctuations around the average
reference value. Using $r_{\mathrm{fit}} = 0.5\lambda_0$ and
$r_{\mathrm{fit}} = 1\lambda_0$ leads to $V = 2.014 \pm 0.034$~km/s and
$V = 2.008 \pm 0.021$~km/s, respectively, so the $V_0 = 2$~km/s
reference value is well recovered with an uncertainty level of 1.5\%
and~1\%. Figures~\ref{fig:C_control}(c) and~(d)
display the spatial error distribution obtained from
\mbox{Equation}~(\ref{eq:error_V}), where the values indicate a similar
error range. Both illustrations demonstrate that an increasing data
range improves accuracy and precision. Another strong feature are the
boundary effects in Figures~\ref{fig:C_control}(a)
and~(b). Recall that we only use observations inside
the square target area, similar to an array deployment. The estimation
error increases towards the boundaries and the affected area appears to
depend on the $r_{\mathrm{fit}}$ length. However, the biasing effects
are not homogeneously distributed along the boundaries but are largest
along the southern edge. This is likely explained by the centered, low
position of the target area in the cavity, together with the excitation
of specific modes in the cavity that sustain a non-isotropic energy
flux in the scattered wave field. This also explains the corresponding
pattern of larger uncertainties in the lower half of 
Figure~\ref{fig:C_control}(c).{\break} The boundary effects emerge in
all our discussed cases but do not affect our conclusions. These
effects can be mitigated by a different configuration or geometry, or
by averaging over more sources located at more diverse locations.

\begin{figure*}[t!]
\includegraphics{fig03}
\caption{\label{fig:C_control}The control experiment. Focal spot
obtained images for (a)~$r_{\mathrm{fit}} = 0.5\lambda_0$ and
(b)~$r_{\mathrm{fit}} = 1\lambda_0$ for a homogeneous medium with wave
speed $V = 2$~km/s. The reference wavelength $\lambda_0$ is indicated
by the black line and the $r_{\mathrm{fit}}$ range used for the
nonlinear regression is the radius of the black circle. The indicated
wave speeds $V$ are the domain average values and the standard
deviation quantify the fluctuations across the domain. Panels~(c)
and~(d) illustrate the spatial distribution of wave speed errors
obtained with Equation~(\ref{eq:error_V}) for $r_{\mathrm{fit}} =
0.5\lambda_0$ and $r_{\mathrm{fit}} = 1\lambda_0$, respectively.}
{\vspace*{-.5pt}}
\end{figure*}

\subsection{Resolution of the interface between{\hfill\break}
half-spaces} \label{sec:res_halfspaces}

Figure~\ref{fig:C_HS} displays the input wave speed distribution
(Figure~\ref{fig:C_HS}a) and the focal spot-based images of two
homogeneous half-space media for the data ranges $r_{\mathrm{fit}} =
0.5\lambda_0$ (Figure~\ref{fig:C_HS}c) and $r_{\mathrm{fit}} =
1\lambda_0$ (Figure~\ref{fig:C_HS}d). As for the control experiment,
the distributions show small residual fluctuations around the well
resolved input values (the distributions in the left half in panels~(c)
and~(d) yield $V = 2.200 \pm 0.099$~km/s and $V = 2.196 \pm
0.095$~km/s, respectively), and again larger edge effects along the
lower boundary.

\begin{figure*}
\includegraphics{fig04}
\caption{\label{fig:C_HS}Results for the two half-spaces. (a)~The input
model. (b)~Reference (blue) and focal spot based velocity profiles
obtained with two data ranges (red, black) that are averaged along the
$N$-axis. Panels~(c) and~(d) show focal spot based images of the
velocity. The reference wavelength $\lambda_0$ is indicated by the
black line and the $r_{\mathrm{fit}}$ range used for the nonlinear
regression is the diameter of the black circle.}
\end{figure*}

More interesting are the profiles across the domain shown in
Figure~\ref{fig:C_HS}(b). These profiles are averages along the
$N$-axis. The imaged distributions across the interface do not follow
the blue input step function but are spread out. To quantify the
resolution we estimate the width $\Delta L$ of the transition. This is
the distance between the samples where the amplitude of the estimated
profile equals the 5\% and 95\% values of the 0.2~km/s velocity jump,
i.e., when the values are 2.01~km/s and 2.19~km/s.
Figure~\ref{fig:LatRes_HS}(a) enlarges the area around the interface
located at $3.75\lambda_0$. It shows four profiles obtained with
$r_{\mathrm{fit}}$ values ranging from $0.5\lambda_0$ to $2\lambda_0$
in $0.5\lambda_0$ intervals. It confirms the previous observation that
the overall velocity in each half-space is correctly estimated for
every data range. As for the results in Figures~\ref{fig:C_control}
and~\ref{fig:C_HS} we can perhaps discern a weak trend to overestimate
the reference values. We attribute this to the imperfect simulation
configuration since similar effects are not observed in numerical
time-reversal experiments, even in the presence of noise
\citep{Giammarinaro2023}. The width $\Delta L$ is shown in
Figure~\ref{fig:LatRes_HS}(a) as dotted lines at the bottom, with the
vertical dashed lines being the 5\% and 95\% boundaries of the
transition zone. The calculated widths are compiled in
Figure~\ref{fig:LatRes_HS}(b). For $r_{\mathrm{fit}} = 0.5\lambda_0$,
the transition zone has a width of $\Delta L = 0.4\lambda_0$. For
$r_{\mathrm{fit}} = 1\lambda_0$ we obtain $\Delta L = 1.2\lambda_0$,
and for $r_{\mathrm{fit}} = 2\lambda_0$, it is $\Delta L =
2.2\lambda_0$. The spreading effect quantified as the transition width
between the two media can hence well be approximated to scale with the
data range, $\Delta L \approx r_{\mathrm{fit}}$.
Figure~\ref{fig:LatRes_HS}(c) shows the average $r_{\mathrm{fit}}$
dependent $\epsilon_V$ profiles. The lack of features at the position
of the interface indicates that the 10\% velocity contrast does not
influence the reconstruction. The profiles indicate the tendency of an
overall better reconstruction in the right slower domain. This is
likely controlled by the proportionally smaller data range in the left
faster domain, which is a consequence of the constant
$r_{\mathrm{fit}}$ value tied to the reference wavelength in the right
domain. As expected, this trend weakens with increasing
$r_{\mathrm{fit}}$. Figures~\ref{fig:LatRes_HS}(a)
and~(c) together show that an increasing
$r_{\mathrm{fit}}$ leads to a generally increased precision as
$\epsilon_V$ decreases, but to a loss in accuracy around the position
of the \mbox{interface}.

\begin{figure*}
\includegraphics{fig05}
\caption{\label{fig:LatRes_HS}The dependence of the transition width
$\Delta L$ on the data range $r_{\mathrm{fit}}$. (a)~Averaged input and
imaged wave speed profiles using four different data ranges. Averaging
is performed along the $N$-axis. The vertical dashed lines indicate the
location where the amplitudes of the empirical profiles equal the 5\%
and 95\% values of the 0.2~km/s velocity jump. The dotted lines
indicate the estimates of the transition width $\Delta L$. (b)~The
transition width $\Delta L$ as a function of the data distance. (c)~The
averaged error estimate~$\epsilon_V$. Colors in all panels correspond
to the same $r_{\mathrm{fit}}$ values.}
\end{figure*}

\subsection{Resolution and contrast of circular{\hfill\break}
inclusions} \label{sec:res_inclusions}

We now examine the resolution of pairs of $1\lambda_0$-wide circular
inclusions separated by variable distances. The velocity in the
inclusions increases and decreases by 0.5~km/s with respect to the
background reference $V_0 = 2$~km/s. Nonlinear regression of the focal
spot data are performed using data ranges $r_{\mathrm{fit}} =
0.5\lambda_0$ and $r_{\mathrm{fit}} = 1\lambda_0$.
Figure~\ref{fig:C_ResInc} collects in the left two columns the input
wave speed distributions, and the focal spot obtained images for
$r_{\mathrm{fit}} = 0.5\lambda_0$ and $r_{\mathrm{fit}} = 1\lambda_0$.
As a general observation, every inclusion is visible on the images for
both data ranges, and for the three different distances separating the
inclusions. This highlights the effectiveness of the Rayleigh wave
focal spot method for ``discovery mode'' applications \citep{Tsai2023}.
However, the contrast---the difference in ``velocity
amplitude''---depends on the data range. Using $r_{\mathrm{fit}} =
1\lambda_0$ increases the area of the biased velocity estimates. The
quantitative aspect of the method decreases with an increase of the
data range, which appears as an averaging effect, consistent with the
observations across the interface in the half-space case. Panels in the
right two columns in Figure~\ref{fig:C_ResInc} show the data from
Figures~\ref{fig:C_ResInc}(a) to~(f) normalized by
the input wave speed maps (Figures~\ref{fig:C_ResInc}a,
b). The overall neutral or white backgrounds
illustrate the well resolved reference value. The fact that circular
features are visible in Figures~\ref{fig:C_ResInc}(i)
to~(l) demonstrates that the velocity estimates
around the edges are \mbox{biased}. The width of the halo or $\Delta L$
scales again with the data distance, the reconstruction benefits from
smaller $r_{\mathrm{fit}}$ values{\break} (Figures~\ref{fig:C_ResInc}i,
j).

\begin{figure*}
{\vspace*{-.2pc}}
\includegraphics{fig06}
{\vspace*{-3pt}}
\caption{\label{fig:C_ResInc}Results of the circular inclusion
resolution test. The left two columns panels~(a--f) show absolute
values, the right two columns panels~(g--l) show scaled values.
(a,b)~Input models of the three circular inclusion pairs separated by
$1\lambda_0$, $0.5\lambda_0$, and $0.25\lambda_0$. The inclusion
diameter is $1\lambda_0$. Focal spot images obtained with (c,d)
$r_{\mathrm{fit}}=0.5\lambda_0$ and (e,f) 
$r_{\mathrm{fit}}=1\lambda_0$. (g,h)~Normalized input wave speed map.
(i,j)~Images in panels~(c,d) scaled by the corresponding input
distributions in panels~(a,b). (k,l)~Images in panels~(e,f) scaled by
the corresponding input distributions in panels~(a,b). In the lower two
rows, the reference wavelength $\lambda_0$ is indicated by the black
line and the $r_{\mathrm{fit}}$ range used for the nonlinear regression
is the diameter of the black circle.}
{\vspace*{-.3pc}}
\end{figure*}

These interpretations are supported by wave speed profiles across the
inclusions (Figure~\ref{fig:C_LineResInc}), which demonstrate again
that the focal spot image quality, i.e., resolution and contrast,
depends on the data range. Inclusions separated by $0.25\lambda_0$
(Figures~\ref{fig:C_LineResInc}a, b) are well
imaged for small $r_{\mathrm{fit}} = 0.5\lambda_0$. The wave speed
value away from an interface is well approximated. This applies to the
stiff and the compliant inclusions. For $r_{\mathrm{fit}} = 1\lambda_0$
the contrast cannot be recovered, which is linked to the data range
dependent transition zone width that is here similar to the inclusion
diameter. The same mechanism applies to the results obtained with
$r_{\mathrm{fit}} = 1.5\lambda_0$. Importantly, two inclusions can
always be discriminated when the separation distance is
$0.25\lambda_0$, which is yet more obvious from
Figure~\ref{fig:C_ResInc}. For a separation distance of $0.5\lambda_0$
this quality of the imaging approach increases, and inclusions are
completely discriminated for a distance of $1\lambda_0$. Taken
together, the data range dependent sensitivity of the resolution mostly
affects the contrast estimate, but less so the ability to discriminate
two objects. The power to discriminate inclusions separated by only
$0.25\lambda_0$ suggests that the property of super-resolution applies,
similar to focal spot medical imaging results obtained with ultrasound
wavelengths \citep{Zemzemi2020}. Focal spot imaging resolution thus
benefits from high station density across short data ranges. The wave
speed estimates are, however, sensitive to noisy data at short
distances. In turn, longer data distances reduce fluctuations, which
leads to a trade-off between accuracy or contrast and resolution or
discrimination \citep{Giammarinaro2023}.\looseness=-1
\pagebreak

\begin{figure*}
\includegraphics{fig07}
\caption{\label{fig:C_LineResInc}Cross-sections through the circular
inclusions. Profiles of input (blue) and estimated wave speeds for
$r_{\mathrm{fit}} = 0.5\lambda_0$ (orange), $r_{\mathrm{fit}} =
1\lambda_0$ (green), and $r_{\mathrm{fit}} = 1.5\lambda_0$ (red) for
each pair of the circular inclusions. Results are obtained along the
$E$-axis passing through the center of the inclusions. (a,b)~Inclusions
separated by $1\lambda_0$. (c,d)~Inclusions separated by
$0.5\lambda_0$. (e,f)~Inclusions separated by $0.25\lambda_0$.}
{\vspace*{-2.5pt}}
\end{figure*}

\subsection{Imaging random media} \label{sec:res_random}
{\vspace*{-.3pc}}

The last experiment explores focal spot imaging results of
heterogeneous media that are characterized by a von~Karman spectral
density probability function (Equations~(\ref{eq:Xi}), (\ref{eq:P})).
Figure~\ref{fig:V_estimation} shows input wave speed distributions and
the images obtained with data ranges $r_{\mathrm{fit}} = 0.5\lambda_0$
and $r_{\mathrm{fit}} = 1\lambda_0$ together with the input-normalized
images and the formal uncertainties obtained with
Equation~(\ref{eq:error_V}). For each case, we calculate a coefficient
of correlation $R$ between the input distribution and the image. As in
the \mbox{previous} experimental configurations, the overall pattern of
the velocity variations is well retrieved. Peak $R$ correlation values
are obtained at zero-lag, which implies no systematic phase change and
hence the effectiveness of focal spot ``inference-mode'' imaging
\citep{Tsai2023}. The smallest coefficient of correlation is $R = 0.66$
for the small scale, high contrast wave speed variation Medium~1 ($a =
0.5\lambda_0$, $\kappa = 0.1$) and the longer data range
$r_{\mathrm{fit}} = 1\lambda_0$. The best estimation is obtained for
Medium~9 ($a = 2\lambda_0$, $\kappa = 0.6$) imaged with
$r_{\mathrm{fit}} = 0.5\lambda_0$, which corresponds to analysis using
shorter ranges than the correlation length~$a$. This leads to a high
coefficient of correlation $R = 0.97$. Increasing $r_{\mathrm{fit}}$ to
$1\lambda_0$ yields $R = 0.95$. This again confirms our previous
results, the images become smoother with increasing $r_{\mathrm{fit}}$,
which is synonymous with a loss in contrast and hence details. As an
example, in row~6, the green 45-degree trending feature around position
$x = 3\lambda$, $y = 2.5\lambda$ illustrates a narrow low-velocity zone
that is overestimated by the averaging, longer $r_{\mathrm{fit}} =
1\lambda_0$. The other way around, the smoother the input distribution
for larger $a$ and $\kappa$ (Medium~3, 6, 9), the better is the
estimate. The corresponding normalized results in rows~3 and~6
illustrate the imaging quality in a complementary style. Predominantly
neutral distributions indicate overall good reconstructions, the
red-green pattern amplitude range indicates the loss in accuracy, and
the color-speckle size is related to $r_{\mathrm{fit}}$. As an example,
in row~6, the green 45-degree trending feature around position $x =
3\lambda$, $y = 2.5\lambda$ illustrates a \mbox{narrow} low-velocity
zone that is overestimated by the averaging, longer $r_{\mathrm{fit}} =
1\lambda_0$. The spatial variations of the formal uncertainty estimates
in row~4 indicate a medium and wave field dependence. As for the
control experiment (Figure~\ref{fig:C_control}c) the generally larger
error in the lower region is governed by the configuration
(Figure~\ref{fig:cavity}). The resolution of the comparatively large
error spot at $x \approx 6\lambda$, $y \approx 3\lambda$ in the short
correlation length media~1 to~6 highlight the advantage of the local
error estimation. Again as for the control experiment, the larger data
range in row~7 reduces the uncertainty associated with
fluctuations{\break} significantly.

\begin{figure*}
\vspace*{2.5pt}
\includegraphics{fig08}
\caption{\label{fig:V_estimation}Results from the random media case.
The top row shows the heterogeneous velocity distributions that are
synthesized using Equations~(\ref{eq:Xi}) and~(\ref{eq:P}). Values of
the corresponding tuning parameters are given in
Table~\ref{tab:RandomParameters}. Rows~2--4 and rows~5--7 show focal
spot based results obtained with $r_{\mathrm{fit}} = 0.5\lambda_0$ and
$r_{\mathrm{fit}} = 1\lambda_0$, respectively. For each case, the top
rows~2 and~5 show the estimated wave speed distributions, the center
rows~3 and~6 show the ratio between the estimated and the input wave
speed, and the bottom rows~4 and~7 show the local wave speed error
estimates. For each panel, the vertical axis and the horizontal axis
correspond respectively to the North and East axis, scaled by the mean
input wavelength. The $R$ coefficient is the correlation between the
input and the estimated map.}
\vspace*{2pt}
\end{figure*}

We estimate the relative wave speed change $\xi$ from the wave speed
images by inverting Equation~(\ref{eq:Cmap})
{\begin{equation}
\xi(\mathbf{x})=V/V_0-1,
\end{equation}}\unskip
with $V_0 = 2$~km/s. We compile histograms of the $\xi$ distributions
to compare properties of the reference input distributions and of the
obtained images. Figure~\ref{fig:histograms} collects the $\xi$
histograms from the input and from the estimates for $r_{\mathrm{fit}}
= 0.5\lambda_0$ and $r_{\mathrm{fit}} = 1\lambda_0$. The similarity
between reference and image is better for $r_{\mathrm{fit}} =
0.5\lambda_0$ than $r_{\mathrm{fit}} = 1\lambda_0$. The best result is
obtained for the smoothest Medium~9. The green histograms obtained with
$r_{\mathrm{fit}} = 1\lambda_0$ are more narrow and have a higher peak
value around small $\xi$ values compared to the orange
$r_{\mathrm{fit}} = 0.5\lambda_0$ results. This indicates again the
low-pass filter property of larger data ranges observed in the previous
experiments.{\vspace*{-2pt}}

\begin{figure*}
\includegraphics{fig09}
\caption{\label{fig:histograms}Histograms of the relative wave speed
change $\xi$ for the nine media imaged in
Figure~\ref{fig:V_estimation}. Blue data correspond to the input
reference values, and orange and green data correspond to images
obtained with $r_{\mathrm{fit}} = 0.5\lambda_0$ and $r_{\mathrm{fit}} =
1\lambda_0$, respectively.}
\end{figure*}

\section{Discussion} \label{sec:discussion}

Resolution can mean different things in different imaging contexts,
including the number or density of measuring points, the spectral
sensitivity of an imaging device or method, the ability to detect or
discriminate features and to accurately estimate their properties,
contrast in brightness or color, or phase \mbox{fidelity}
\citep[e.g.,][]{Smith2013}. We use numerical simulations of
two-dimensional acoustic wave propagation in a cavity
(Figures~\ref{fig:cavity}, \ref{fig:Example_proc}) to investigate the
lateral resolution of the focal spot imaging technique for a fixed
acquisition system with a constant number of grid points. The increase
in depth resolution for such a compact dense array configuration
compared to measurements made on traveling seismic Rayleigh surface
waves is established by \citet{Giammarinaro2023}. We implement four
test cases that together allow us to study the lateral resolution power
of Rayleigh wave focal spots reconstructed from vertical--vertical
component noise correlation data. Lateral resolution is discussed as
the ability to resolve a step function in the material properties, to
discriminate and characterize two closely spaced objects, and to
measure position, amplitude, and phase of random distributions.

In a first homogeneous control experiment
(Figures~\ref{fig:Example_proc}, \ref{fig:C_control}) the reference
wave speed is estimated with an error below 2\% on average, which
includes, however, areas of larger bias associated with edge effects.
This average error is larger compared to the focal spot results based
on numerical time-reversal experiments, where noise-free synthetics
lead to errors in the 0.01\% range for vertical--vertical component
records and data ranges of $r_{\mathrm{fit}} = 0.5\lambda_0$ and
$r_{\mathrm{fit}} = 1\lambda_0$, and where anisotropic surface wave
incidence results in biases in the 0.1\% range
\citep{Giammarinaro2023}. The different error levels are associated
with the different methods used to synthesize Green's functions and
focal spots. In time-reversal experiments, the wave field and hence the
ballistic wave correlations are fully controlled by the mirror
properties. More mirror elements lead to better refocusing results.
Controlled lab experiments can stack over different space realizations.
In seismic data applications an improved Green's function and
refocusing reconstruction is achieved by time averaging to better
conform with the decorrelated noise source assumption. Here, the focal
spots are retrieved from cross-correlation of the reverberating cavity
wave field, where the quality of the Green's function is controlled by
the ability to excite and average a sufficiently large number of modes
in the cavity \citep{Draeger1999}. Hence we stack over different
realizations using different source positions. This increases the
number of independent modes to enhance the narrow-band refocusing,
which is equivalent to using more time reversal mirror elements in a
time-reversal experiment. Our approach converges towards the
theoretical focal spot, but it remains sensitive to details of the
implementation such as the source positions of the relatively few
realizations. This explains the observed fluctuations in the focal spot
reconstructions, which indicate the imperfect Green's function
synthesis, even for the homogeneous case. This means that our imperfect
cavity results are comparable to focal spots obtained from noisy data. 

\looseness=-1
The second half-spaces experiment (Figures~\ref{fig:C_HS},
\ref{fig:LatRes_HS}) shows the feasibility to resolve the contrast
between two media that have a 10\% difference in wave speed. However,
the interface is not perfectly resolved. Whereas the wave speed is
correctly estimated away from the interface, the finite data or fitting
range creates an averaging or low-pass filter effect that depends on
$r_{\mathrm{fit}}$. Tests with different $r_{\mathrm{fit}}$ values
indicate that the transition width $\Delta L$ scales with good
approximation linearly with the data range, $\Delta L \approx
r_{\mathrm{fit}}$.

The third circular inclusion configuration (Figures~\ref{fig:C_ResInc},
\ref{fig:C_LineResInc}) shows the possibility to discriminate two
separate objects separated by $0.25\lambda_0$, even for data ranges of
one wavelength. However, the data range dependent low-pass filter
properties lead to averaged amplitude values at edges. We have not
observed situations where the focal spot method yields biased phase
properties, so the inclusion positions are accurate. This suggests that
the super-resolution property demonstrated with passive elastography in
soft tissues \citep{Zemzemi2020} also applies to 2D Rayleigh surface
wave propagation. Again this means that for good data quality and high
station density the method has the potential to meet the formal
criterion of super-resolution, i.e., the sensitivity at small scales is
sufficient to discriminate objects or features that are separated by
distances that are much shorter than the wavelength.

The fourth test case considers random media
(Figures~\ref{fig:V_estimation}, \ref{fig:histograms},
Table~\ref{tab:RandomParameters}) which show the highest similarity to
distributions of Earth material properties. The quality of the focal
spot reconstruction as quantified by the correlation coefficient $R$
between reference and image depends on the roughness or smoothness of
the distribution in relation to wavelength and data range. $R$ is small
when the distributions are rough, have small-scale fluctuations
compared to the probing wavelength, and the data range is large. The
reconstruction is almost perfect when the distributions are smooth,
have large-scale fluctuations, and the data range is small. These
conclusions are further corroborated comparing histogram properties of
the velocity variation parameter $\xi$ (Equation~(\ref{eq:Xi}),
Figure~\ref{fig:histograms}), which again demonstrate the low-pass
filtering effect of large $r_{\mathrm{fit}}$ values. Thus, positions
are well estimated but the amplitudes diverge for small-scale
heterogeneities. The best estimates are obtained for variations on
scales larger than $r_{\mathrm{fit}}$. Deblurring methods can
potentially be developed to further improve the lateral resolution of
the method considering that the focal spot Bessel function shape acts
as a spatial convolution filter. If resources permit, an imaging
campaign could be optimized by first detecting target features with low
contrast before then a sensor re-configuration helps to improve the
quantitative estimate by increasing the signal-to-noise ratio through
network densification. 

The present study employs two-dimensional acoustic simulations. These
scalar simulations yield the same Green's function as for
vertical--vertical component Rayleigh wave propagation. Our main
conclusions thus hold for seismic Rayleigh wave imaging, including the
advantageous estimation of the local error that can help improve the
uncertainty management of shear wave inversions. However, this set-up
does not allow to study the biasing effect of interfering body wave
energy. \citet{Giammarinaro2023} shows that the error on the
vertical--vertical component increases in the presence of $P$~waves,
but this can be compensated for by increasing the data range. Together
with our observations here this implies that an increase in
$r_{\mathrm{fit}}$ improves the results if refocusing $P$~wave energy
distorts the surface wave focal spot, but that this remedy negatively
affects lateral resolution power. This trade-off situation would
benefit from efficient focal spot filters. Alternatively, Rayleigh wave
phase speed estimates obtained from radial-vertical component data are
much less sensitive to $P$~waves \citep{Giammarinaro2023}, which offers
independent constraints for the improvement of vertical--vertical
results. The application of radial--radial and transversal--transversal
component focal spots requires first an efficient separation of the
combined Rayleigh and Love wave energy \citep{Haney2014}. Spatial
autocorrelations of surface wave strain and rotational data can further
provide additional constraints on the local velocity structure
\citep{Nakahara2021, Nakahara2022}.

\section{Conclusion} \label{sec:conclusion}

We investigate the lateral resolution power of the Rayleigh surface
wave focal spot imaging method using two-dimensional numerical
experiments of reverberating wave fields in a cavity. Most importantly
the resolution depends on the data range. This means that focal spot
imaging exhibits super-resolution properties provided the data quality
supports sub-wavelength data ranges. Longer data ranges still allow
imaging of small-scale features at super-resolution albeit with a loss
in contrast. Seismic Rayleigh wave focal spot imaging shows convincing
resolution properties that make it suitable for a wide range of imaging
applications ranging from feature detection to accurate wave speed
estimates. There are hence no fundamental disadvantages compared to
established passive surface wave tomography methods. Here as there, the
station configuration can be tuned to support image quality and
properties for different goals, and in both cases the data quality or
signal-to-noise ratio ultimately has the largest impact on the
resolution.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\section*{Open research/data availability}

No observational data were used in this study.

\section*{Dedication}

The manuscript was written with contributions from all authors. All
authors have given approval to the final version of the manuscript.

\section*{Funding}
This work was sponsored by a research grant from the Academy of Finland
(decision number 322421).

\section*{Acknowledgments}

We thank Julien de~Rosny, Stefan Catheline, Michel Campillo, L\'eonard
Seydoux, Laurent Stehly, Alexander Meaney, and Markus Juvonen for
helpful discussions. We thank the editor N.~Shapiro and two anonymous
reviewers for comments that helped to improve the manuscript.

\CDRGrant[Academy of Finland]{322421}

\back{}

\bibliographystyle{crgeos}
\bibliography{crgeos20221200}
\refinput{crgeos20221200-reference.tex}

\end{document}