[a]Yi Chen

# First measurement of anti-kT jet spectra and jet substructure using the
archived ALEPH $e^{+}e^{-}$ data at 91.2 GeV

Austin Baty  Dennis Perepelitsa  Christopher McGinn  Jesse Thaler  Marcello
Maggi  Paoti Chang  Tzu-An Sheng  Yang-Ting Chien  Yen-jie Lee

###### Abstract

We present the first anti-kT jet spectrum and substructure measurements using
the archived ALEPH $e^{+}e^{-}$ data taken in 1994 at a center of mass energy
of $\sqrt{s}=91.2$ GeV. Jets are reconstructed with the anti-kT algorithm with
a resolution parameter of 0.4. It is the cleanest test of jets and QCD without
the complication of hadronic initial states. The fixed center-of-mass energy
also allows the first direct test of pQCD calculation. We present both the
inclusive jet energy spectrum and the leading dijet energy spectra, together
with a number of substructure observables. They are compared to predictions
from PYTHIA6, PYTHIA8, Sherpa, HERWIG, VINCIA, and PYQUEN. None of the models
fully reproduce the data. The data are also compared to two perturbative QCD
calculations at NLO and with NLL’+R resummation. The results can also serve as
reference measurements to compare to results from hadronic colliders. Future
directions, including testing jet clustering algorithms designed for future
electron-ion collider experiments, will also be discussed.

##  1 Introduction

Experimental signatures similar to those typically attributed to the Quark-
gluon plasma have been observed in progressively smaller systems, for example
the ridge-like enhancement in two-particle long range correlation in high
multiplicity proton-proton collisions at the LHCs [1]. Jet-related
observables, however, do not show any signs of modification so far in smaller
systems. It is therefore of interest for jet measurements in the smallest
$e^{+}e^{-}$ collision systems. Previous jet measurements in LEP (see e.g.
[2]) were limited to earlier generations of jet-finding algorithms and
couldn’t be directly compared with recent results from hadronic collisions at
the LHC and RHIC (e.g. [3, 4, 5, 6]), where the anti-kT [7] algorithm is
commonly used. In the current work [8], we performed the first measurement of
energy spectrum and substructure of anti-kT jet using the archived data from
the ALEPH experiment.

##  2 Jet reconstruction and calibration

Detector-level jets are clustered [9] from the energy-flow objects with the
anti-kT algorithm with a distance parameter of 0.4, using the spherical-
coordinate variant of the algorithm suitable for the $e^{+}e^{-}$ collision
system with energy $E$ instead of transverse momentum $k_{\mathrm{T}}$ and
opening angle $\theta$ instead of $\Delta
R\equiv\sqrt{\Delta\eta^{2}+\Delta\phi^{2}}$. Due to the small inactive area
around the beam pipe, in order to avoid jets that overlap with it, only jets
within the acceptance region $0.2\pi<\theta_{\text{jet}}<0.8\pi$ are
considered. A hadronic event selection [10] is applied to reject events with
only electromagnetic interactions.

Jets are calibrated with a multi-step procedure. As the first step, they are
corrected to truth-level jet energies using the archived simulated Pythia6
[11] sample in different bins of jet direction $\theta_{\text{jet}}$.
Additional data-simulation difference are derived. Through leading dijet
energies the $e^{+}$-going and the $e^{-}$-going sides of the detector are
intercalibrated. The absolute scale is then derived using the event-wide
multijet mass, comparing data to simulation, excluding tails which may contain
additional effects. The magnitude of the data-simulation difference goes up to
1%.

The energy resolution of jets are evaluated first in simulation. The
difference in resolution between data and simulation are measured using the
leading dijet energy balance, using the third-leading jet as a handle for
systematics on the method. The measured relative resolution difference is up
to 5% between simulation and data.

For the jet substructure considered in this work, the SoftDrop [12, 13]
grooming algorithm is used with parameter $(z_{\text{cut}}=0.1,\beta=0.0)$.
The algorithm cleans soft particles scattered at large angles and identifies
two subjets, which can be used to define observables. A few observables are
measured: $z_{G}$, defined as the energy fraction carried by the smaller of
the two subjets; $R_{G}$, the opening angle between the two subjets;
$M_{G}/E$, the invariant mass divided by the two subjets normalized by the
energy of the jet; and $M/E$, the invariant mass of the whole jet before the
grooming algorithm normalilzed by the energy of the jet. Results with
additional parameter settings and different algorithms are left for future
studies.

##  3 Analysis

The results are unfolded to account for detector smearing effects. One example
transfer matrix is shown in Fig. 1 for measurement of jet energy spectrum. The
difference in resolution between data and simulation is incorporated. Due to
the migration in energy between truth-level and detector-level, for
substructure measurements the unfolding is done in two dimensions in
conjunction with the jet energy.

Figure 1: Transfer matrix for inclusive jet energy measurement.

In order to better quantify the in-cone energy, spectrum of global leading
dijet energy is measured. A correction is derived from simulation to account
for cases where some of the energies of the jet is missing due to overlap with
the beam direction.

Sources of systematic uncertainties include jet energy scale and resolution,
contribution from combinatorial jets not associated with a parton from the
hard process, unfolding, and modeling. For the energy spectra measurements,
the dominant source is jet energy scale and resolution, while for the
substructure measurements the dominant is from the modeling.

##  4 Results

The unfolded energy spectrum is shown in Fig. 2. In the left panel the
measured energy spectrum is compared with Pythia6 [11], Pythia8 [14], Herwig
7.2.2 [15] and Sherpa [16] generators. The peak above 40 GeV comes from the
dominant $e^{+}e^{-}\rightarrow q\bar{q}$ process, while the rise at lower
energy are from multijet topologies or energies escaped from the leading jets.
The generators generally describe the peak region well, but some discrepancies
are observed at low energy. The data is also compared with NLL’-resummed
calculation [17], and a general good agreement is observed. Parton level NLO
spectrum is also shown.

Figure 2: Measured jet energy spectrum compared with Monte-Carlo simulations
(left) and NLL’-resummed calculations and parton-level NLO spectrum (right).

In Fig. 3 the measured leading dijet energy is shown. Even though some trends
are seen in the ratio plots, the larger systematic uncertainty prevents a
precise statement. In this case the dominant uncertainty comes from the
modeling uncertainty originating from the fact that only one set of detector-
level simulation is available to use for studies.

Figure 3: Measured leading dijet jet energy compared with Monte-Carlo
simulations (left) and NLL’-resummed calculations (right).

The groomed jet energy balance $z_{G}$ and radius $R_{G}$ are shown in Fig. 4
for jets with $E>40$ GeV. Jets that do not result in two valid subjets through
the SoftDrop algorithm are shown in the bin below 0. The $z_{G}$ distribution
for high energy jets is similar to what has been observed in other collision
systems. There is a deviation up to 10% between the data and simulations,
consistent with the level of discrepancy observed in proton-proton collisions
at the LHC. The bulk of the $R_{G}$ distribution is well-described by
simulations, but there are some discrepancies at the tail of the $R_{G}$
distribution.

Figure 4: Measured groomed jet energy sharing $z_{G}$ (left) and $R_{G}$
(right) for inclusive jets with jet $E>40$ GeV. The data is compared with
predictions from various generators.

An example of groomed jet mass normalized by jet energy ($M_{G}/E$) is shown
in Fig. 5. In the left panel jets with small energy is shown, while in the
right panel we show jets with $E>40$ GeV. The general behavior is distinct
between low and high energies, indicating that the particles are more spread
out for low energy jets, whereas the higher energy jets behave more similar to
the high $p_{T}$ jets seen in hadron colliders. The general shape is
nevertheless captured by simulations. Similar to the observation in $R_{G}$,
the tail of the distributions are less well captured by the simulations.

Figure 5: Measured groomed jet mass normalized by jet energy for jets with
$10<E<15$ GeV (left) and $E>40$ GeV (right) for inclusive jets. The data is
compared with predictions from various generators.

The full list of results, including results from other jet energy ranges, full
jet mass $M/E$ and the sum of leading dijet energy, can be found in ref [8].

##  5 Summary

We presented the first measurement of the energy and substructure of anti-kT
jets in $e^{+}e^{-}$ collisions using the archived ALEPH data. The data is
generally described by popular event generators, however, some deviations are
present for low jet energy and substructure observables. The result can
provide input to QCD theory calculations and event generators. The calibrated
jets also serve as an excellent testing ground for new algorithms in a clean
environment without hadronic initial states.

## References

  * [1] Vardan Khachatryan et al.  Observation of Long-Range Near-Side Angular Correlations in Proton-Proton Collisions at the LHC.  JHEP, 09:091, 2010.
  * [2] D. Buskulic et al.  Quark and gluon jet properties in symmetric three jet events.  Phys. Lett., B384:353–364, 1996.
  * [3] A. Abdesselam et al.  Boosted Objects: A Probe of Beyond the Standard Model Physics.  Eur. Phys. J. C, 71:1661, 2011.
  * [4] Gavin P. Salam.  Towards Jetography.  Eur. Phys. J. C, 67:637–686, 2010.
  * [5] A. Altheimer et al.  Jet Substructure at the Tevatron and LHC: New results, new tools, new benchmarks.  J. Phys. G, 39:063001, 2012.
  * [6] Andrew J. Larkoski, Ian Moult, and Benjamin Nachman.  Jet Substructure at the Large Hadron Collider: A Review of Recent Advances in Theory and Machine Learning.  Phys. Rept., 841:1–63, 2020.
  * [7] Matteo Cacciari, Gavin P. Salam, and Gregory Soyez.  The anti-$k_{t}$ jet clustering algorithm.  JHEP, 04:063, 2008.
  * [8] Yi Chen et al.  Jet energy spectrum and substructure in $e^{+}e^{-}$ collisions at 91.2 GeV with ALEPH archived data.  JHEP, 2022(6):008, 2022.
  * [9] Matteo Cacciari, Gavin P. Salam, and Gregory Soyez.  FastJet User Manual.  Eur. Phys. J. C, 72:1896, 2012.
  * [10] Anthony Badea, Austin Baty, Paoti Chang, Gian Michele Innocenti, Marcello Maggi, Christopher Mcginn, Michael Peters, Tzu-An Sheng, Jesse Thaler, and Yen-Jie Lee.  Measurements of two-particle correlations in $e^{+}e^{-}$ collisions at 91 GeV with ALEPH archived data.  Phys. Rev. Lett., 123(21):212002, 2019.
  * [11] Torbjörn Sjöstrand, Patrik Eden, Christer Friberg, Leif Lonnblad, Gabriela Miu, Stephen Mrenna, and Emanuel Norrbin.  High-energy physics event generation with PYTHIA 6.1.  Comput. Phys. Commun., 135:238–259, 2001.
  * [12] Andrew J. Larkoski, Simone Marzani, Gregory Soyez, and Jesse Thaler.  Soft Drop.  JHEP, 05:146, 2014.
  * [13] Mrinal Dasgupta, Alessandro Fregoso, Simone Marzani, and Gavin P. Salam.  Towards an understanding of jet substructure.  JHEP, 09:029, 2013.
  * [14] Torbjörn Sjöstrand, Stefan Ask, Jesper R. Christiansen, Richard Corke, Nishita Desai, Philip Ilten, Stephen Mrenna, Stefan Prestel, Christine O. Rasmussen, and Peter Z. Skands.  An introduction to PYTHIA 8.2.  Comput. Phys. Commun., 191:159–177, 2015.
  * [15] Johannes Bellm et al.  Herwig 7.0/Herwig++ 3.0 release note.  Eur. Phys. J. C, 76(4):196, 2016.
  * [16] T. Gleisberg, Stefan. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert, and J. Winter.  Event generation with SHERPA 1.1.  JHEP, 02:007, 2009.
  * [17] Duff Neill, Felix Ringer, and Nobuo Sato.  Leading jets and energy loss.  JHEP, 07:041, 2021.


# Note on stability of an abstract coupled hyperbolic-parabolic system:
singular case

Kaïs Ammari  LR Analysis and Control of PDEs, LR 22ES03, Department of
Mathematics, Faculty of Sciences of Monastir, University of Monastir, Tunisia
<EMAIL_ADDRESS>,  Farhat Shel  LR Analysis and Control of PDEs, LR
22ES03, Department of Mathematics, Faculty of Sciences of Monastir, University
of Monastir, Tunisia<EMAIL_ADDRESS>and  Zhuangyi Liu  Department of
Mathematics and Statistics, University of Minnesota, Duluth, MN 55812-3000,
United States<EMAIL_ADDRESS>
###### Abstract.

In this paper we try to complete the stability analysis for an abstract system
of coupled hyperbolic and parabolic equations

$\left\\{\begin{array}[]{lll}\displaystyle u_{tt}+Au-A^{\alpha}w=0,\\\
w_{t}+A^{\alpha}u_{t}+A^{\beta}w=0,\\\
u(0)=u_{0},u_{t}(0)=u_{1},w(0)=w_{0},\end{array}\right.$

where $A$ is a self-adjoint, positive definite operator on a complex Hilbert
space $H$, and $(\alpha,\beta)\in[0,1]\times[0,1]$, which is considered in
[1], and after, in [4]. Our contribution is to identify a fine scale of
polynomial stability of the solution in the region
$S_{3}:=\left\\{(\alpha,\beta)\in[0,1]\times[0,1];\,\beta<2\alpha-1\right\\}$
taking into account the presence of a singularity at zero.

###### Key words and phrases:

Hyperbolic-parabolic system, stability

###### 2010 Mathematics Subject Classification:

35B65, 35K90, 47D03

###### Contents

  1. 1 Introduction
  2. 2 Stabilization
  3. 3 Application

##  1\. Introduction

In this paper, we study the stability of the following system :

(1.1) $u_{tt}+Au-A^{\alpha}w=0,$ (1.2) $w_{t}+A^{\alpha}u_{t}+A^{\beta}w=0,$
(1.3) $u(0)=u_{0},u_{t}(0)=u_{1},w(0)=w_{0},$

where $A$ is a self-adjoint, positive definite operator on a complex Hilbert
space $H$, and $(\alpha,\beta)\in
S_{3}=\left\\{(a,b)\in[0,1]\times[0,1];\,b<2a-1\right\\}$.

By denoting $U=(u,u_{t},w)^{T},U_{0}=(u_{0},u_{1},w_{0})^{T}$, system
(1.1)-(1.3) can be written as an abstract linear evolution equation on the
space $\mathcal{H}=\mathcal{D}(A^{\frac{1}{2}})\times H\times H$,

(1.4)
$\left\\{\begin{array}[]{ll}\displaystyle\frac{dU}{dt}(t)=\mathcal{A}_{\alpha,\beta}U(t),\,t\geq
0,\\\ U(0)=U_{0},\end{array}\right.$

where the operator
$\mathcal{A}_{\alpha,\beta}:\mathcal{D}(\mathcal{A}_{\alpha,\beta})\subset\mathcal{H}\rightarrow\mathcal{H}$
is defined by

$\mathcal{A}_{\alpha,\beta}=\left(\begin{array}[]{ccll}0&I&0\\\
-A&0&A^{\alpha}\\\ 0&-A^{\alpha}&-A^{\beta}\end{array}\right),$

with the domain

$\mathcal{D}(\mathcal{A}_{\alpha,\beta})=\mathcal{D}(A)\times\mathcal{D}(A^{\alpha})\times\mathcal{D}(A^{\alpha}),$

Firstly, we give a detailed review about the well-posedness of the problem
(1.4). The operator $\mathcal{A}_{\alpha,\beta}$ is densely defined and
dissipative, we will prove that its closure generates a $C_{0}$-semigroup of
contractions. Using a Lumer-Phillips theorem [6], it suffices to prove that
the adjoint operator $\mathcal{A}^{*}_{\alpha,\beta}$ is also dissipative.

The operator $\mathcal{A}^{*}_{\alpha,\beta}$ is a closed extension of the
operator

$\mathcal{M}_{\alpha,\beta}=\left(\begin{array}[]{ccll}0&-I&0\\\
A&0&-A^{\alpha}\\\ 0&A^{\alpha}&-A^{\beta}\end{array}\right),$

with domain

$\mathcal{D}(\mathcal{M}_{\alpha,\beta})=\mathcal{D}(\mathcal{A}_{\alpha,\beta})=\mathcal{D}(A)\times\mathcal{D}(A^{\alpha})\times\mathcal{D}(A^{\alpha}).$

The operator $\mathcal{M}_{\alpha,\beta}$ is densely defined and dissipative,
then it is closable and its $\overline{\mathcal{M}}_{\alpha,\beta}$ is also
dissipative. To conclude, it suffices to prove that
$\mathcal{A}^{*}_{\alpha,\beta}=\overline{\mathcal{M}}_{\alpha,\beta}$ . For
this we use the following lemma [2, 1].

######  Lemma 1.1.

[2] We consider on the Hilbert spaces $G$ and $H_{1}$ the operators

$\displaystyle\mathcal{A}:\mathcal{D}(\mathcal{A})\subset G\rightarrow
G,\;\;\;\;\;B:\mathcal{D}(B)\subset H_{1}\rightarrow G,$ $\displaystyle
B^{*}:\mathcal{D}(B^{*})\subset G\rightarrow
H_{1},\;\;\;\;\;\mathcal{C}:\mathcal{D}(\mathcal{C})\subset H_{1}\rightarrow
H_{1}.$

and we consider the operator matrix $\mathcal{M}$ on $G\times H_{1}$ defined
by

$\mathcal{M}:=\left(\begin{array}[]{ccll}\mathcal{A}&B\\\
-B^{*}&\mathcal{C}\end{array}\right),\;\;\;\mathcal{D}(\mathcal{M}):=\left(\mathcal{D}(\mathcal{A})\cap\mathcal{D}(B^{*})\right)\times\left(\mathcal{D}(B)\cap\mathcal{D}(\mathcal{C})\right).$

Assume that $\mathcal{C}$ is boundedly invertible, that
$\mathcal{B}\in\mathcal{L}\left(\mathcal{D}(\mathcal{C}),G\right)$ and that
$\mathcal{C}^{-1}\mathcal{B}^{*}$ extends to a bounded linear operator (we
denote its closure with the same symbol). Then $\mathcal{M}$ is closed if and
only if $\mathcal{A}+\mathcal{B}\mathcal{C}^{-1}\mathcal{B}^{*}$ is closed and
one has

$\displaystyle\overline{\mathcal{M}}=\left(\begin{array}[]{ccll}I&\mathcal{B}\mathcal{C}^{-1}\\\
0&I\end{array}\right)\left(\begin{array}[]{ccll}\overline{\mathcal{A}+\mathcal{B}\mathcal{C}^{-1}\mathcal{B}^{*}}&0\\\
0&\mathcal{C}\end{array}\right)\left(\begin{array}[]{ccll}I&0\\\
-\mathcal{C}^{-1}\mathcal{B}^{*}&I\end{array}\right)$
$\displaystyle\mathcal{D}\left(\overline{\mathcal{M}}\right)=\left\\{\left(\begin{array}[]{cll}u\\\
v\end{array}\right)\in\mathcal{D}(\overline{\mathcal{A}+\mathcal{B}\mathcal{C}^{-1}\mathcal{B}^{*}})\times
H_{1},\;-\mathcal{C}^{-1}\mathcal{B}^{*}u+v\in\mathcal{D}(\mathcal{C})\right\\}.$

By taking $G=\mathcal{D}(A^{1/2})$, $H_{1}=H\times H$, $\mathcal{A}=0$,
$\mathcal{B}=(-I\;0):\mathcal{D}(\mathcal{B})=\mathcal{D}(A^{\alpha})\times\mathcal{D}(A^{\alpha})\subset
H\times H\rightarrow\mathcal{D}(A^{1/2})$,
$\mathcal{C}=\left(\begin{array}[]{ccll}0&-A^{\alpha}\\\
A^{\alpha}&-A^{\beta}\end{array}\right)$, with
$\mathcal{D}(\mathcal{C})=\mathcal{D}(A^{\alpha})\times\mathcal{D}(A^{\alpha})$,
it appears that $\mathcal{B}^{*}=\left(\begin{array}[]{cll}-A\\\
0\end{array}\right)$,
$\mathcal{C}^{-1}\mathcal{B}^{*}=\left(\begin{array}[]{cll}A^{1+\beta-2\alpha}\\\
A^{1-\alpha}\end{array}\right)$,
$\mathcal{A}+\mathcal{B}\mathcal{C}^{-1}\mathcal{B}^{*}=-A^{\beta+1-2\alpha}$
(with domain $\mathcal{D}(A)$), and that
$\mathcal{D}\left(\overline{\mathcal{A}+\mathcal{B}\mathcal{C}^{-1}\mathcal{B}^{*}}\right)=\mathcal{D}(A^{1/2})$.
Furthermore, the operators $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$
satisfies the hypothesis in Lemma 1.1, then

$\mathcal{D}\left(\overline{\mathcal{M}}_{\alpha,\beta}\right)=\left\\{(u,v,\theta)^{T}\in\mathcal{H},\;-A^{1+\beta-2\alpha}u+v\in\mathcal{D}(A^{\alpha}),\;-A^{1-\alpha}u+\theta\in\mathcal{D}(A^{\alpha})\right\\}.$

But a direct calculation, gives

$\mathcal{D}\left(\mathcal{A}^{*}_{\alpha,\beta}\right)\subset\left\\{(u,v,\theta)^{T}\in\mathcal{H},\;-A^{1+\beta-2\alpha}u+v\in\mathcal{D}(A^{\alpha}),\;-A^{1-\alpha}u+\theta\in\mathcal{D}(A^{\alpha})\right\\}.$

Hence $\mathcal{A}^{*}_{\alpha,\beta}=\overline{\mathcal{M}}_{\alpha,\beta}$
and $\overline{\mathcal{A}}_{\alpha,\beta}$ generates a $C_{0}$-semigroup of
contractions.

Now, concerning the stability of the semigroup
$e^{t\mathcal{A}_{\alpha,\beta}}$, recall that Ammar-Khodja et al. [1], are
firstly proved that for every $\alpha,\beta\geq 0$,
$e^{t\mathcal{A}_{\alpha,\beta}}$ is exponentially stable if and only if
$\max(1-2\alpha,2\alpha-1)<\beta<2\alpha.$ Later on, a stability analysis has
been performed by J. Hao and Z. Liu in [5] for
$(\alpha,\beta)\in[0,1]\times[0,1]$. They divided the unit square
$[0,1]\times[0,1]$ into four regions $S$, $S_{1}$, $S_{2}$, $S_{3}$ where

$\displaystyle S$ $\displaystyle:=$
$\displaystyle\left\\{(\alpha,\beta)\in[0,1]\times[0,1];\,\max(1-2\alpha,2\alpha-1)\leq\beta\leq
2\alpha\right\\},$ $\displaystyle S_{1}$ $\displaystyle:=$
$\displaystyle\left\\{(\alpha,\beta)\in[0,1]\times[0,1];\;0<\beta-2\alpha,\;\alpha\geq
0,\;\frac{1}{2}\leq\beta\leq 1\right\\},$ $\displaystyle S_{2}$
$\displaystyle:=$
$\displaystyle\left\\{(a,b)\in[0,1]\times[0,1];\;\beta<1-2\alpha,\;\alpha\geq
0,\;0\leq\beta\leq\frac{1}{2}\right\\}$

and $S_{3}$ as defined below. We summarize the main results in the following
theorem (see also [4]).

######  Theorem 1.2.

The semigroup $e^{t\mathcal{A}_{\alpha,\beta}}$ has the following stability
properties:

(i) In $S$, it is exponentially stable;

(ii) In $S_{1}$, it is polynomially stable of order
$\frac{1}{2(\beta-2\alpha)}$;

(iii) In $S_{2}$, it is polynomially stable of order
$\frac{1}{2-2(\beta+2\alpha)}$;

(iv) In $S_{3}$, it is not asymptotically stable.

To justify the non asymptotic stability in region $S_{3}$, they proved that
$0\in\sigma(\mathcal{A}_{\alpha,\beta})$, where
$\sigma(\mathcal{A}_{\alpha,\beta})$ is the spectrum of
$\mathcal{A}_{\alpha,\beta}$ . Moreover, it can be shown that
$\sigma(\mathcal{A}_{\alpha,\beta})\cap\mathbf{i}\mathbb{R}=\left\\{0\right\\}$
in the region $S_{3}$.

The main result of this paper then concerns the precise asymptotic behaviour
of the solutions of (1.1)-(1.3) for $(\alpha,\beta)$ in the region $S_{3}$,
with initial condition in a special subspace of
$\mathcal{D}(\mathcal{A}_{\alpha,\beta})$. Precisely, we will estimate
$\|e^{t\mathcal{A}_{\alpha,\beta}}\mathcal{A}_{\alpha,\beta}\left(I-\mathcal{A}_{\alpha,\beta}\right)^{-1}\|$
when $t\longrightarrow\infty$. Our technique is a special frequency and
spectral analysis of the corresponding operator.

##  2\. Stabilization

We will justify that the resolvent has only a singularity (at zero) on the
imaginary axis (Lemma 2.3 below), and that
$\left(\mathbf{i}\lambda-\mathcal{A}_{\alpha,\beta}\right)^{-1}$ is bounded
outside a neighborhood of zero in $\mathbb{R}$ (Lemma 2.4 below). then we
apply a result due to Batty, Chill and Tomilov ([3, Theorem 7.6 ]) which
relate the decay of
$\|e^{t\mathcal{A}_{\alpha,\beta}}\mathcal{A}_{\alpha,\beta}\left(I-\mathcal{A}_{\alpha,\beta}\right)^{-1}\|$
to the growth of
$\left(\mathbf{i}\lambda-\mathcal{A}_{\alpha,\beta}\right)^{-1}$ near zero.

So, first we recall the corresponding result

######  Theorem 2.1.

([3], Theorem 7.6.) Let $(T(t))_{t\geq 0}$ be a bounded $C_{0}$-semigroup on a
Hilbert space $X$, with generator $\mathcal{A}$. assume that
$\sigma(\mathcal{A})\cap\mathbf{i}\mathbb{R}=\left\\{0\right\\}$, and let
$\gamma\geq 1$. The following are equivalent:

(i)
$\|\left(\mathbf{i}s-\mathcal{A}\right)^{-1}\|=\left\\{\begin{array}[]{ll}\displaystyle
O\,(|s|^{-\gamma}),\;\;\,s\rightarrow 0,\\\ \displaystyle
O\,(1),\;\;\;\;\;\;\;\;\,|s|\rightarrow\infty,\end{array}\right.$

(ii)
$\|T(t)A^{\gamma}(I-A)^{-\gamma}\|=O\,\left(\frac{1}{t}\right),\;\;\;\;t\rightarrow\infty,$

(iii)
$\|T(t)A(I-A)^{-1}\|=O\,\left(\frac{1}{t^{1/\gamma}}\right),\;\;\;\;t\rightarrow\infty$.

Hence, our main result is the following

######  Theorem 2.2.

We have the following decay in the region $S_{3}$

(2.1)
$\|e^{t\mathcal{A}_{\alpha,\beta}}\mathcal{A}_{\alpha,\beta}\left(I-\mathcal{A}_{\alpha,\beta}\right)^{-1}\|=O\,\left(\frac{1}{t}\right),\;\;\;\;t\rightarrow\infty.$

In particular, the decay of
$e^{t\mathcal{A}_{\alpha,\beta}}\mathcal{A}_{\alpha,\beta}x$ to zero is
uniform with respect to $x\in\mathcal{D}(\mathcal{A}_{\alpha,\beta})$.

It follows also that, for every $z\in Ran(\mathcal{A}_{\alpha,\beta})$ we have

$\|e^{t\mathcal{A}_{\alpha,\beta}}z\|=O\,\left(\frac{1}{t}\right),\;\;\;\;t\rightarrow\infty.$

###### Proof.

In view of Theorem 2.1, the proof is a direct consequence of the following
three lemmas.

######  Lemma 2.3.

In the region $S_{3}$, the operator $\mathcal{A}_{\alpha,\beta}$ satisfies

$\sigma(\mathcal{A}_{\alpha,\beta})\cap\mathbf{i}\mathbb{R}=\left\\{0\right\\}.$

###### Proof.

It has been proved in [5], at the third section that
$\left\\{0\right\\}\subset\sigma(\mathcal{A}_{\alpha,\beta})\cap\mathbf{i}\mathbb{R}$.

Conversely, to show that there is no nonzero spectrum point on the imaginary
axis, we use a contradiction argument. In fact, let $\lambda\in\mathbb{R}$,
$\lambda\neq 0$ such that
$\mathbf{i}\lambda\in\sigma(\mathcal{A}_{\alpha,\beta})$. Then, there exists a
sequence $(U_{n})\subset\mathcal{D}(\mathcal{A}_{\alpha,\beta})$, with
$\|U_{n}\|=1$ for all $n$, such that

(2.2) $\lim_{n\rightarrow\infty}\|\left(\mathbf{i}\lambda
I-\mathcal{A}_{\alpha,\beta}\right)U_{n}\|=0$

or there exists a sequence
$(U_{n})\subset\mathcal{D}(\mathcal{M}_{\alpha,\beta})=\mathcal{D}(\mathcal{A}_{\alpha,\beta})$,
with $\|U_{n}\|=1$ for all $n$, such that

(2.3) $\lim_{n\rightarrow\infty}\|\left(\mathbf{i}\lambda
I+\mathcal{M}_{\alpha,\beta}\right)U_{n}\|=0.$

Setting $U_{n}=(u_{n},v_{n},\theta_{n})$, then (2.2) is equivalent to

(2.4) $\displaystyle\mathbf{i}\lambda
A^{1/2}u_{n}-A^{1/2}v_{n}=o(1),\;\;\;\text{in}\;H,$ (2.5)
$\displaystyle\mathbf{i}\lambda
v_{n}+Au_{n}-A^{\alpha}\theta_{n}=o(1),\;\;\;\text{in}\;H,$ (2.6)
$\displaystyle\mathbf{i}\lambda\theta_{n}+A^{\alpha}v_{n}+A^{\beta}\theta_{n}=o(1),\;\;\;\text{in}\;H.$

First, since

$Re\left(\left\langle\left(\mathbf{i}\lambda
I-\mathcal{A}_{\alpha,\beta}\right)U_{n},U_{n}\right\rangle_{\mathcal{H}}\right)=\|A^{\beta/2}\theta_{n}\|^{2}$

we obtain

(2.7) $\lim_{n\rightarrow\infty}\|A^{\beta/2}\theta_{n}\|=0,$

and in particular

(2.8) $\lim_{n\rightarrow\infty}\|\theta_{n}\|=0.$

Second, taking inner product of (2.4) with $\frac{1}{\lambda}A^{1/2}u_{n}$,
(2.5) with $\frac{1}{\lambda}v_{n}$ and (2.6) with
$\frac{1}{\lambda}\theta_{n}$, taking into account (2.7) and (2.8), we get

(2.9)
$\displaystyle\mathbf{i}\|A^{1/2}u_{n}\|^{2}-\frac{1}{\lambda}\left\langle
v_{n},Au_{n}\right\rangle=o(1),$ (2.10)
$\displaystyle\mathbf{i}\|v_{n}\|^{2}+\frac{1}{\lambda}\left\langle
Au_{n},v_{n}\right\rangle-\frac{1}{\lambda}\left\langle
A^{\alpha}\theta_{n},v_{n}\right\rangle=o(1),$ (2.11)
$\displaystyle\frac{1}{\lambda}\left\langle
A^{\alpha}v_{n},\theta_{n}\right\rangle=o(1).$

Then, by combining (2.9)-(2.11) and using that $\|U_{n}\|=1$, it yields

(2.12) $\|A^{1/2}u_{n}\|=\frac{1}{2}+o(1),\;\;\;\;\|v_{n}\|=\frac{1}{2}+o(1).$

Now, since in the region $S_{3}$, $1-\alpha<\frac{1}{2}$ and
$\frac{3}{2}-2\alpha+\beta<\frac{1}{2}$, then using the boundedness of
$\|A^{1/2}u_{n}\|$ we get, by interpolation,

(2.13) $\|A^{1-\alpha}u_{n}\|=O(1),\;\;\;\;\|A^{1-2\alpha+\beta}u_{n}\|=O(1).$

Then, replacing $\mathbf{i}v_{n}$ by $\lambda u_{n}$ in (2.5), due to (2.4)
and interpolation, next taking the inner product of the obtained equation with
$\frac{1}{\lambda}A^{1-2\alpha+\beta}u_{n}$ to get

(2.14)
$-\lambda\|A^{\frac{1}{2}-\alpha+\frac{\beta}{2}}u_{n}\|^{2}+\frac{1}{\lambda}\|A^{1-\alpha+\frac{\beta}{2}}u_{n}\|^{2}-\frac{1}{\lambda}\left\langle
A^{\alpha}\theta_{n},A^{1-2\alpha+\beta}u_{n}\right\rangle=o(1).$

Taking the inner product of (2.6) with $\frac{1}{\lambda}A^{1-\alpha}u_{n}$,
we get

(2.15)
$\mathbf{i}\left\langle\theta_{n},A^{1-\alpha}u_{n}\right\rangle+\frac{1}{\lambda}\left\langle
A^{1/2}v_{n},A^{1/2}u_{n}\right\rangle+\frac{1}{\lambda}\left\langle
A^{\alpha}\theta_{n},A^{1-2\alpha+\beta}u_{n}\right\rangle=o(1).$

By (2.8) and (2.13), the first term in (2.14) converge to zero. Moreover,
using (2.4), we can replace $\frac{1}{\lambda}A^{1/2}v_{n}$ in the second term
in (2.15) by $\mathbf{i}A^{1/2}u_{n}$. Consequently, the sum of (2.14) and
(2.15) yields

(2.16)
$-\lambda\|A^{\frac{1}{2}-\alpha+\frac{\beta}{2}}u_{n}\|^{2}+\frac{1}{\lambda}\|A^{1-\alpha+\frac{\beta}{2}}u_{n}\|^{2}+\mathbf{i}\|A^{\frac{1}{2}}u_{n}\|^{2}=o(1).$

Hence, $\|A^{\frac{1}{2}}u_{n}\|^{2}=o(1)$, which contradict the first
estimate in (2.12).

The same approach applied to (2.3) leads to the same conclusion without any
difficulty. ∎

######  Lemma 2.4.

In the region $S_{3}$,

$\|\left(\mathbf{i}s-\mathcal{A}_{\alpha,\beta}\right)^{-1}\|=O\,(1),\;\;\,s\rightarrow\infty.$

###### Proof.

It is a direct consequence of Theorem 2.3 in [5], since
$\limsup\limits_{\lambda\in\mathbb{R},\,\lambda\rightarrow\infty}|\lambda|^{\frac{\beta}{\alpha}}\|\left(\mathbf{i}\lambda-\mathcal{A}_{\alpha,\beta}\right)^{-1}\|<\infty$.
∎

######  Lemma 2.5.

In the region $S_{3}$,

$\|\left(\mathbf{i}s-\mathcal{A}_{\alpha,\beta}\right)^{-1}\|=O\,(|s|^{-1}),\;\;\,s\rightarrow
0.$

###### Proof.

By contradiction, suppose that $\limsup\limits_{s\in\mathbb{R},\,s\rightarrow
0}\|s\left(\mathbf{i}sI-\mathcal{A}_{\alpha,\beta}\right)^{-1}\|=\infty$.

Putting $s=\frac{1}{w}$, this is equivalent to
$\limsup\limits_{w\in\mathbb{R},\,|w|\rightarrow\infty}\|w^{-1}\left(\mathbf{i}w^{-1}I-\mathcal{A}_{\alpha,\beta}\right)^{-1}\|=\infty$.
Then, there exists a sequence $(w_{n})$ of real numbers with
$|w_{n}|\rightarrow\infty$, and a sequence
$(U_{n})\subset\mathcal{D}(\mathcal{A}_{\alpha,\beta})$, with $\|U_{n}\|=1$
for all $n$ such that

(2.17)
$\lim_{n\rightarrow\infty}|w_{n}|\|\left(\mathbf{i}w_{n}^{-1}I-\mathcal{A}_{\alpha,\beta}\right)U_{n}\|=0.$

Setting $U_{n}=(u_{n},v_{n},\theta_{n})$, then (2.17) is rewretten explicitely
as follows

(2.18)
$\displaystyle\mathbf{i}w_{n}^{-1}|w_{n}|A^{1/2}u_{n}-|w_{n}|A^{1/2}v_{n}=o(1),\;\;\;\text{in}\;H,$
(2.19)
$\displaystyle\mathbf{i}w_{n}^{-1}|w_{n}|v_{n}+|w_{n}|Au_{n}-|w_{n}|A^{\alpha}\theta_{n}=o(1),\;\;\;\text{in}\;H,$
(2.20)
$\displaystyle\mathbf{i}w_{n}^{-1}|w_{n}|\theta_{n}+|w_{n}|A^{\alpha}v_{n}+|w_{n}|A^{\beta}\theta_{n}=o(1),\;\;\;\text{in}\;H.$

Since
$Re\left(\left\langle|w_{n}|\left(\mathbf{i}w_{n}^{-1}I-\mathcal{A}_{\alpha,\beta}\right)U_{n},U_{n}\right\rangle_{\mathcal{H}}\right)=-|w_{n}|\|A^{\beta/2}\theta_{n}\|^{2}$,
it yields

(2.21) $\lim_{n\rightarrow\infty}|w_{n}|^{1/2}\|A^{\beta/2}\theta_{n}\|=0,$

and in particular

(2.22) $\lim_{n\rightarrow\infty}\|\theta_{n}\|=0.$

Then, applying $|w_{n}|^{-1}A^{-\alpha}$ to (2.20), we get, (by taking into
account (2.22)),

(2.23) $\lim_{n\rightarrow\infty}\|v_{n}\|=0.$

Now, taking inner product of (2.18) with $A^{1/2}u_{n}$, (2.19) with $v_{n}$
and (2.20) with $\theta_{n}$, taking into account (2.21), (2.22) and (2.23),
we have

(2.24)
$\displaystyle\mathbf{i}w_{n}^{-1}|w_{n}|\|A^{1/2}u_{n}\|^{2}-|w_{n}|\left\langle
v_{n},Au_{n}\right\rangle=o(1),$ (2.25) $\displaystyle|w_{n}|\left\langle
Au_{n},v_{n}\right\rangle-|w_{n}|\left\langle
A^{\alpha}\theta_{n},v_{n}\right\rangle=o(1),$ (2.26)
$\displaystyle|w_{n}|\left\langle
A^{\alpha}\theta_{n},v_{n}\right\rangle=o(1).$

By combining (2.24)-(2.26), it yields

(2.27) $\lim_{n\rightarrow\infty}\|A^{1/2}u_{n}\|=0.$

The promised contradiction follows from (2.22), (2.23) and (2.27). Thus, the
proof of Lemma 2.5 is completed. ∎

Which ends the proof of theorem. ∎

######  Remark 2.6.

In $S_{3}$, the semigroup $e^{t\mathcal{A}_{\alpha,\beta}}$ is of Gevrey class
$\delta>\frac{\alpha}{\beta}$ and in particular, it is infinitely
differentiable [5]. Thus, for $z\in Ran(\mathcal{A}_{\alpha,\beta})$, we not
only has polynomial stability, but also have instantaneous smoothness.

##  3\. Application

As application of Theorem 2.2, we consider the following thermoplate system:

$\left\\{\begin{array}[]{lll}\displaystyle
u_{tt}+\Delta^{2}u-(-\Delta)^{\alpha}w=0,\Omega\times(0,+\infty),\\\
w_{t}+(-\Delta)^{\alpha}u_{t}-\Delta w=0,\,\Omega\times(0,+\infty),\\\
u=\Delta u=w=0,\,\Gamma\times(0,+\infty),\\\
u(x,0)=u_{0}(x),u_{t}(x,0)=u_{1}(x),w(x,0)=w_{0}(x),\,\Omega,\end{array}\right.$

where $\Omega$ be a bounded domain in $\mathbb{R}^{n}$ with smooth boundary
$\Gamma$, $\beta=\frac{1}{2}$ and $3/4<\alpha\leq 1.$

Here, $H=L^{2}(\Omega),A=\Delta^{2}$ avec $\mathcal{D}(A)=\left\\{u\in
H^{2}(\Omega)\cap H^{1}_{0}(\Omega);\,\Delta u=0\;\text{on}\;\Gamma\right\\}$.

Then, according to Theorem 2.2 the corresponding semigroup satisfies the
estimate (2.1) for all $\alpha\in(3/4,1].$

## References

  * [1] F. Ammar-Khodja, A. Bader and A. Benabdallah, Dynamic stabilization of systems via decoupling techniques, ESAIM Control Optim. Calc. Var., 4 (1999), 577–593.
  * [2] F. V. Atkinson, H. Langer, R. Mennicken and A. A. Shkalikov, The essential spectrum of some matrix operators. Math. Nachr., 167 (1994), 5–20.
  * [3] C. J. K. Batty, R. Chill and Y. Tomilov, Fine scales of decays of operator semigroups. J. Eur. Math. Soc., 18 (2016), 853–929. stability of linear dynamical systems in Hilbert space, Ann. Differential Equations, 1 (1985), 43–56.
  * [4] J. Hao and Z. Liu, Stability of an abstract system of coupled hyperbolic and parabolic equations, Z. Angew. Math. Phys., 64 (2013), 1145–1159.
  * [5] J. Hao, Z. Liu and J. Yong, Regularity analysis for an abstract system of coupled hyperbolic and parabolic equations, Journal of Differential Equations, 259 (2015), 4763-4798.
  * [6] G. Lumer and R. S. Phillips, Dissipative operators in a Banach space, Pacific Journal of Mathematics, 11 (1961), 679-698.


# Spectral Gap-Based Seismic Survey Design

Oscar López1, Rajiv Kumar2, Nick Moldoveanu3 and Felix J. Herrmann4
1Harbor Branch Oceanographic Institute, Florida Atlantic University
2Schlumberger WesternGeco
3Independent Geophysical Consultant
4School of Earth and Atmospheric Sciences, Georgia Institute of Technology

###### Abstract

Seismic imaging in challenging sedimentary basins and reservoirs requires
acquiring, processing, and imaging very large volumes of data (tens of
terabytes). To reduce the cost of acquisition and the time from acquiring the
data to producing a subsurface image, novel acquisition systems based on
compressive sensing, low-rank matrix recovery, and randomized sampling have
been developed and implemented. These approaches allow practitioners to
achieve dense wavefield reconstruction from a substantially reduced number of
field samples.

However, designing acquisition surveys suited for this new sampling paradigm
remains a critical and challenging role in oil, gas, and geothermal
exploration. Typical random designs studied in the low-rank matrix recovery
and compressive sensing literature are difficult to achieve by standard
industry hardware. For practical purposes, a compromise between stochastic and
realizable samples is needed. In this paper, we propose a deterministic and
computationally cheap tool to alleviate randomized acquisition design, prior
to survey deployment and large-scale optimization. We consider universal and
deterministic matrix completion results in the context of seismology, where a
bipartite graph representation of the source-receiver layout allows for the
respective spectral gap to act as a quality metric for wavefield
reconstruction. We provide realistic scenarios to demonstrate the utility of
the spectral gap as a flexible tool that can be incorporated into existing
survey design workflows for successful seismic data acquisition via low-rank
and sparse signal recovery.

###### Index Terms:

Spectral gap, low rank matrix completion, compressive sensing, seismic data,
seismic trace interpolation, bipartite graph, biadjacency matrix, nuclear norm
minimization.

##  I Introduction

Seismic exploration is a complex and lengthy process involving many expensive
components. To ensure high quality imaging in frugal manner, survey designs
are an important initial step that determine the efficiency of data
acquisition before costly in-field deployment [40, 24, 31]. The typical
components for data acquisition survey designs are:

  1. 1.

Definition of seismic survey objectives in terms of reservoir depth, type of
reservoir, required vertical and horizontal imaging resolution.

  2. 2.

Analysis of the existent seismic, well, geological data, near surface and
environmental conditions in the survey area.

  3. 3.

Seismic modeling and imaging of different acquisition scenarios to evaluate
placement of sources and receivers based on resolution, signal-to-noise
requirements, and inherent noise levels.

  4. 4.

Selection of optimum acquisition parameters based on modeling, imaging
results, and the cost of acquiring and processing the data.

In all steps, it is essential for a survey design study to have a realistic
earth model of the area. If geological and geophysical information in the
survey area do not allow an earth model to be built, existing models available
in the industry can be informative when the imaging challenges of the survey
area are comparable. For instance, the SEG Seam model is quite often use for
survey design studies to evaluate the acquisition geometries for sub-salt
exploration [27].

To further complicate the survey design process, recent years have introduced
a new paradigm for data acquisition due to the advent of compressive sensing
(CS) [10]. Low-rank matrix recovery (LRMR) and CS based seismic exploration
has successfully materialized from an academic novelty to an indispensable
industrial tool [17, 29, 23, 44]. As a popular example, seismic data
reconstruction from limited measurements provides an instance of the _matrix
completion_ problem [33, 4, 6, 5]. These approaches allow for simultaneous
acquisition and compression of seismic data, where a substantially reduced
number of field samples allow dense wavefield reconstruction via large-scale
optimization. As a consequence, practitioners may achieve high-fidelity
imaging and post-processing with economical surveys.

While presenting many advantages, LRMR and CS-based acquisition remains a
challenging task with many open questions in addition to the pre-existing
challenges of Nyquist-based survey design. In particular, it is not clear how
to fulfill step 3) above in an practical manner with respect to LRMR and CS.
In the theory of matrix completion, generating appropriate sampling schemes is
typically guaranteed with high probability by selecting entries to observe
uniformly at random from a dense grid [33, 4, 6, 5]. This is difficult to
achieve in seismology since it corresponds to sources and receivers removed
uniformly at random from equipment arrays independently for each survey.

This incompatibility between theory and practice further complicates step 4)
for non-Nyquist-based acquisition design. Such a process would involve
simulating a 5D data tensor (corresponding to a 3D seismic survey) and
performing large-scale reconstruction for all possible options within physical
constraints, e.g., solving the _nuclear norm minimization_ problem [35]. As a
consequence, LRMR and CS-based acquisition design is a complex procedure that
requires computational exploration of numerous criteria to decide upon an
appropriate source-receiver layout.

To mitigate such tedious processes, we propose a computationally cheap metric
to evaluate a given acquisition design. The novel idea comes from the matrix
completion literature [3, 11, 32], where the authors show that matrix
completion is successful if the spectral gap (SG) of the associated set of
observed entries is large. To elaborate, let $\textbf{X}\in\mathbb{C}^{n\times
m}$ represent densely sampled seismic data (in some transform domain) reshaped
as a matrix that exhibits low-rank structure, i.e., X is well approximated
(according to some tolerance) by a rank-$r$ matrix where $r\ll\min(n,m)$ [17,
23]. Relative to this desired densely sampled data, we subsample and let
$\textbf{M}\in\\{0,1\\}^{n\times m}$ be a binary matrix whose non-zero entries
specify the sampled matrix entries. The SG is defined as the difference
between M’s first and second singular values
$\sigma_{1}(\textbf{M}),\sigma_{2}(\textbf{M})$. The main results in [3, 32]
state that one can recover an incoherent (i.e., “not spiky”) low-rank matrix
via _nuclear norm minimization_ if the number of samples and the SG are large
enough. Intuitively, the SG can be seen as a qualitative measure informing us
how well we can expect a given sampling scheme to perform matrix completion.

In this work, we propose the SG as a deterministic measure to differentiate
appropriate surveys at an initial design stage prior to computationally
expensive simulations and in-field acquisition. This paper argues that the SG
provides a simple but powerful tool that can help severely lessen the amount
of numerical experiments required to evaluate pre-survey field designs (e.g.,
simulations via data models and large-scale optimization). Our work aims to
demonstrate that the SG can be incorporated into survey planning schemes as a
numerically cheap and flexible step to reduce the overall design complexity.

The relationship between the SG and quality of reconstruction is demonstrated
in Figure 1. Given $\textbf{M}\in\\{0,1\\}^{n\times m}$ specifying the sampled
data entries, we constitute its SG by the ratio
$\frac{\sigma_{2}(\textbf{M})}{\sigma_{1}(\textbf{M})}\in[0,1]$. Henceforth,
we will refer to the term
$\frac{\sigma_{2}(\textbf{M})}{\sigma_{1}(\textbf{M})}$ as the _SG ratio_ and
note that a large SG corresponds to a small SG ratio
($\frac{\sigma_{2}(\textbf{M})}{\sigma_{1}(\textbf{M})}\ll 1$) and vice versa.
Figure 1 illustrates the correlation between M’s SG ratio and the signal-to-
noise ratio (SNR, defined in (3)) of reconstruction via LRMR with observed
matrix entries according to M. In each experiment, with fixed sampling
percentage, only the distribution of the observed entries and corresponding SG
change. Therefore, the plot showcases better quality of reconstruction from
samples with larger SG. The full details of the conducted experiments are
postponed until Section III-A. Further details for generating M from a survey
design are discussed in Section II-A.

(a)  Figure 1: Illustration of the relationship between the quality of
reconstruction and the SG ratio of the sampled matrix entries. The figure
plots the average reconstruction SNR (defined in (3)) vs average SG ratio
$\frac{\sigma_{2}(\textbf{M})}{\sigma_{1}(\textbf{M})}$ of the corresponding
binary sampling matrix M, for 3 percentages of missing receivers. The results
are taken from the average of 100 independent experiments.

Design of optimal acquisition geometries suitable for CS is an ongoing area of
research with several solutions in the literature relevant to our work. In
[28], the authors propose a reconstruction based design procedure where a
greedy method outputs an appropriate sensing matrix. The methodology requires
data interpolation at each step of the optimization problem, thus making it a
computationally demanding process. The work in [1] uses prior information from
the low-frequency spectrum of seismic data to construct a signal and coherent
noise model, followed by the evaluation of a quality metric for acquisition
layout as a function of frequency. The optimal acquisition layout corresponds
to the worst case value of the metric over the band. The proposed method does
not require any synthetic or real data, but the metric involves computing
several condition numbers (i.e, ratio of the smallest and largest singular
values of the sensing matrix). This ratio is numerically expensive for sensing
matrices related to large-scale 5D seismic data. In contrast, our proposed
method requires computing the two largest singular values of a single binary
matrix representing source-receiver layout. This approach is numerically
efficient and does not require underlying data nor a prior model.

We validate the SG metric with several numerical experiments that consider
distinct aspects of the seismic data acquisition and reconstruction process.
In the next section, we begin with some background from the literature of
matrix completion that formally introduces the SG. Therein, Section II-A gives
a detailed description of how the spectral gap can be computed from a survey
design. Section III is dedicated to numerical experiments, where in Section
III-A we empirically demonstrate the relationship between the SG and quality
of reconstruction via LRMR. Section III-B demonstrates how the SG is sensitive
to large gaps of missing data and in Section III-C we utilize the SG to
determine an appropriate density of the output interpolation grid. Section
III-D argues that the SG is also informative in sparsity-based recovery. We
end with concluding remarks and future work in Section IV.

##  II Universal and Deterministic Matrix Completion

In contrast to most results in the matrix completion literature, the work in
[3, 11, 32] aims for universal and deterministic recovery guarantees. The term
universal refers to the set of observed matrix entries,
$\Omega\subset\\{1,2,\cdots,n\\}\times\\{1,2,\cdots,m\\}$, where this subset
is expected to provide an accurate output via nuclear norm minimization for an
entire class of signals uniformly. Such results are of great value to
practitioners since $\Omega$ will function for all matrices of interest,
removing the need to design new sampling schemes for every independent seismic
survey. Furthermore, the attached conditions that ensure the suitability of
$\Omega$ inherently provide a means to quantify the source-receiver layout.

This last observation is the focus of this article, where the spectral gap
serves as a computationally cheap step to distinguish between given
acquisition designs. The proposed metric differentiates sampling patterns in
terms of their suitability for reconstruction via matrix completion. When
incorporated into existing survey design schemes, the SG offers an important
tool that can offset numerically expensive simulation-based acquisition
design.

To elaborate, we begin with some matrix completion background. In particular,
we focus on deterministic and universal matrix completion [3, 11, 32] where
the authors quantify the success of matrix completion via the spectral gap. In
what follows, let $n\geq m$ without loss of generality.

Given a matrix of interest $\textbf{X}\in\mathbb{C}^{n\times m}$ and a subset
of observed entries $\Omega\subset\\{1,2,\cdots,n\\}\times\\{1,2,\cdots,m\\}$,
our collected data is given as
$\textbf{B}=P_{\Omega}(\textbf{X}+\textbf{E})\in\mathbb{C}^{n\times m}$, where
for $(k,\ell)\in\\{1,2,\cdots,n\\}\times\\{1,2,\cdots,m\\}$

$P_{\Omega}(\textbf{X}+\textbf{E})_{k\ell}=\left\\{\begin{array}[]{ll}\textbf{X}_{k\ell}+\textbf{E}_{k\ell}&\mbox{if
}(k,\ell)\in\Omega\\\ 0&\mbox{otherwise }\end{array}\right.$ (1)

and E encompasses the noise in our observations with bounded Frobenius norm

$\|P_{\Omega}(\textbf{E})\|_{F}:=\left(\sum_{(k,\ell)\in\Omega}|\textbf{E}_{k\ell}|^{2}\right)^{1/2}\leq\eta.$

Matrix completion via nuclear norm minimization aims to approximate

$\textbf{X}\approx\textbf{X}^{\sharp}:=\underset{\textbf{Z}\in\mathbb{C}^{n\times
m}}{\text{argmin}}\quad\|\textbf{Z}\|_{*}\ \ \mbox{s.t.}\ \
\|P_{\Omega}(\textbf{Z})-\textbf{B}\|_{F}\leq\eta,$ (2)

where

$\|\textbf{Z}\|_{*}=\sum_{k=1}^{m}\sigma_{k}(\textbf{Z})$

is the _nuclear norm_ and $\sigma_{k}(\textbf{Z})$ denotes the $k$-th largest
singular value of Z. This method has been extensively studied [33, 4, 6, 5].
In the case $\eta=0$, standard results in the literature show that X can be
exactly recovered via (2) with high probability under conditions of the
general form:

  * •

rank$(\textbf{X})=r\leq\min\\{n,m\\}$.

  * •

X is incoherent, e.g., satisfies the strong incoherence property (see [33, 4,
6, 5, 3]). Intuitively, this assumption ensures that the energy (i.e.,
Frobenius norm) of the data matrix is not concentrated on a small subset of
entries.

  * •

$\Omega$ is generated uniformly at random with $|\Omega|\sim nr\log^{2}(n)$,
where $|\Omega|$ indicates the number of observed entries. Furthermore, each
matrix to be recovered requires an independent $\Omega$ to be generated.

The first condition is crucial since equation (2) with $|\Omega|<nm$ is an
underdetermined problem and the low-rank constraint aims to compensate for the
limited number of observations. This validates the use of the nuclear norm as
the objective function, in order to output $\textbf{X}^{\sharp}$ that can be
written as a sum of few rank 1 matrices (i.e., low-rank structure). The second
assumption is a necessary condition for sampling schemes that are oblivious of
the matrix structure, but can be avoided if one samples according to prior
knowledge of the matrix’s leverage scores [7].

However, the third condition is restrictive and does not hold in many
applications. For example, in seismic acquisition the allowed structure of
$\Omega$ is limited by measurement hardware and physical constraints involved
in data collection. Achieving uniform random sampling independently for each
survey would require a significant amount of novel sensing equipment. Instead,
matrix completion results that consider achievable sampling schemes would be
most informative and applicable.

To this end, we now consider results involving universal and deterministic
matrix completion [3, 11, 32]. Let 1 be the all-ones $n\times m$ matrix and
$\textbf{M}:=P_{\Omega}(\textbf{1})\in\\{0,1\\}^{n\times m}$. We will refer to
M as the sampling mask and note that it is a binary matrix, which specifies
the locations of sampled matrix entries. The main result in [3] shows that if
the spectral gap (SG) of M is large enough (i.e.,
$\sigma_{2}(\textbf{M})\ll\sigma_{1}(\textbf{M})$) then one can recover any
low-rank incoherent matrix via (2). Specifically, the result reads as follows:

######  Theorem II.1 (Theorem 4.2 in [3]).

Let $\textbf{X}\in\mathbb{C}^{n\times m}$ be a rank-$r$ matrix satisfying the
strong incoherence property with parameter $\mu$ (see [5] and Claim 5.1 in
[3]). Let $\Omega$ be generated such that the top singular vectors of M are
the all 1’s vectors and

$|\Omega|\geq
36\frac{\sigma_{2}(\textbf{M})^{2}}{\sigma_{1}(\textbf{M})}\mu^{2}nr^{2}.$

Then X is the unique optimum of (2) with $\eta=0$.

In the result, the multiplicative factor 36 is an absolute constant (not
necessarily optimal) obtained by the authors, independent of the matrix
dimensions. This factor bears no particular significance and can likely be
improved to further reduce the number of samples needed for universal matrix
completion. The incoherence parameter $\mu$ quantifies how evenly spread the
information is throughout the data matrix and ensures that we sample
accordingly. The condition on the singular vectors of M requires $\Omega$ to
contain the same number of samples in each row and column. This constraint is
quite restrictive, but also required for the main results in [11, 32].
Generating sampling schemes with this structure is not straightforward [19,
20] and complex to translate to field samples. However, in Section III, we
allow ourselves the liberty to violate this assumption. Our numerical
experiments will demonstrate empirically that the relationship between the SG
and the success of matrix completion remains useful even when this condition
does not hold.

In contrast to standard uniform random sampling requirements,
$|\Omega|\sim\mathcal{O}(nr\log^{2}(n))$, Theorem II.1 requires
$|\Omega|\sim\mathcal{O}(nr^{2})$. While this sample complexity is pessimistic
when $r>\log^{2}(n)$, we stress that the result holds deterministically and
universally for all incoherent matrices. In the best case scenario we can
expect
$\frac{\sigma_{2}(\textbf{M})^{2}}{\sigma_{1}(\textbf{M})}\sim\mathcal{O}(1)$
(e.g., Ramanujan graph [13, 3]). More importantly, Theorem II.1 provides a
useful relationship between the success of matrix completion and the ratio
$\frac{\sigma_{2}(\textbf{M})}{\sigma_{1}(\textbf{M})}\in[0,1]$ directly
related to the SG. This ratio appears with similar implications in the related
work [11, 32]. Intuitively, a practitioner choosing amongst several sampling
designs can compute the respective SG ratios and make an educated choice of
which design will output the best quality of matrix reconstruction.

###  II-A LRMR-Based Seismic Data Acquisition

Our work proposes the SG as a useful metric for survey planning that
differentiates sampling schemes at varying stages of existing designs. Keeping
in mind the results from the previous section, our proposed matrix completion-
based seismic data acquisition and reconstruction hinges upon the following
three requisites:

  1. i)

Densely sampled seismic data (in some transform domain) reshaped as a matrix
$\textbf{X}\in\mathbb{C}^{n\times m}$ should exhibit incoherence and low-rank
structure, i.e., X is well approximated (according to some tolerance) by a
rank-$r$ matrix where $r\ll\min(n,m)$.

  2. ii)

A computationally efficient program for large-scale LRMR, e.g., nuclear norm
minimization [35].

  3. iii)

The subset of observed matrix entries
$\Omega\subset\\{1,2,\cdots,n\\}\times\\{1,2,\cdots,m\\}$ should be suitable
for incoherent matrices.

To address i), we note that several reorganizations and transformations of
seismic data tensors into incoherent matrices have been shown to expunge low-
rank structure [17, 23]. Furthermore, for ii) efficient methodologies have
been proposed to solve the matrix completion problem with large matrices [2,
18, 15], including on parallel architectures [34, 45].

This paper focuses on fulfilling iii) in a practical manner using the SG.
While specific circumstances will vary, our general approach consists of
generating a sampling mask M for each acquisition under consideration and
comparing them via the SG. To do so, a practitioner must first determine the
desired output sampling density as in steps 1) and 2) from Section I. This
will implicitly dictate the matrix dimensions $n$ and $m$, which also depend
on the transform domain chosen to satisfy requisite i) above. A given sampling
scheme will sub-sample relative to this desired source-receiver spacing, i.e.,
collect wavefield measurements on a subset of
$\\{1,2,\cdots,n\\}\times\\{1,2,\cdots,m\\}$ and thereby determine $\Omega$.
The sampling mask is then defined as the $n\times m$ matrix with value of 1 on
the entries in $\Omega$ and 0’s otherwise. Finally, the singular values
$\sigma_{1}(\textbf{M}),\sigma_{2}(\textbf{M})$ and the respective gap (or SG
ratio) can be computed.

Comparing the SG ratios of each acquisition design considered provides a
computationally efficient mechanism to select or eliminate sub-sampling
schemes based on their suitability for matrix completion. This procedure can
alleviate the complexity of simulation-based evaluations for survey design,
where the SG can reduce the number of computationally intense experiments
typically required for steps 3) and 4) from Section I.

It is important to notice that in the context of seismology, the SG is
oblivious to the underlying geological complexity. Therefore, the SG tool
proposed here does not consider many important attributes of seismic survey
design related to specific models or geological scenes (e.g., evaluating
azimuth coverage and fold [42, 14]). Such topics are an important component of
standard acquisition design [25], and should be addressed in conjunction. In
this sense, our proposed SG metric should not be considered sufficient for a
complete survey design but rather as a tool that can reduce the numerical
complexity of existing workflows. In Section III, we demonstrate how the SG
can be used to limit the amount of data-dependent evaluations and how it
largely succeeds at evaluating surveys under models of varying geological
complexity.

###  II-B Connections with Graph Theory

Interestingly, the SG is a common tool for analyzing the connectivity of
communication networks [13, 43]. In this context, a larger SG provides a more
robust network. The foundation comes from graph theory and the study of
expander graphs, where a communication network can be seen as a group of
vertices (network nodes) with corresponding edges (signifying communication
between nodes). Considering the respective adjacency matrix (analogous to M in
our context), this matrix’s SG gives a direct measurement of the network’s
efficiency in the flow of information.

A similar interpretation can be presented in our seismology scenario. The
array of sources and receivers can be seen as the vertices with edges that
represent the recording and transmission of pressure waves. In other words,
using the matricization in Section III-A (see also [16, 17]), a sourcex-
receiverx pair (vertex in the $x$-axis) will have an edge with a sourcey-
receivery pair (vertex in the $y$-axis) if the receiverx-receivery pair
recorded the pressure wave due to the corresponding sourcex-sourcey pair (and
vice versa with respect to the reciprocal role of transmission). Thus, in our
case the sampling mask M can be seen as the corresponding biadjacency matrix
of the bipartite graph that captures a sampling scenario with such vertices
and edges. Since the interpolation problem implicitly depends on well
distributed samples, we can intuitively see that a well connected graph (i.e.,
large SG) will be more effective at infilling data between source-receiver
pairs. This interpretation was first elaborated abstractly in [4], where the
authors show that observed entries generating a fully connected graph are
necessary to recover a matrix of any rank.

##  III Numerical Experiments

To demonstrate the utility of the SG for survey design, we numerically address
four specific questions related to seismic data acquisition and
reconstruction:

  * •

Section III-A: What is the relationship between the SG and the quality of
LRMR-based reconstruction?

  * •

Section III-B: Is the SG sensitive to large gaps of missing data?

  * •

Section III-C: Can the SG identify appropriate output interpolation grid
density?

  * •

Section III-D: Is the SG informative for sparsity-based reconstruction?

We design numerical experiments to provide insight for these questions. Some
of our experiments are abstract in the sense that we do not consider explicit
seismic data, and instead focus on the source-receiver distribution. This
stresses that the SG can be applied prior to data collection and without
simulated data. However, most industrial acquisition designs require
simulation on data models [41, 39]. To demonstrate the utility of the SG in
such workflows, we also provide experiments with validation on geologically
diverse synthetic marine datasets (BG 3D, SEG/EAGE Overthrust and the SEAM
models).

These experiments aim to demonstrate that even though the SG is oblivious to
the underlying data, the reconstruction quality predicted by this tool largely
pans out for models of varying geological complexity. While these examples do
not prove the proposed method will always work, they indicate that the
spectral gap is a good metric for survey suitability that can be incorporated
into standard simulation-based design methodologies. Furthermore, we note that
our experiments mainly focus on seismology in marine settings. Land surveys
require distinct considerations and separate experiments are needed to fully
explore the role of the SG in such acquisitions. However, due to its
detachment from the underlying data, we believe that the experiments here will
be informative for the applicability of the SG in land acquisition geometries.

In what follows, we use the signal-to-noise ratio (SNR) as a metric of matrix
reconstruction quality defined as

$-20\log_{10}\left(\frac{\|\textbf{X}^{\sharp}-\textbf{X}\|_{F}}{\|\textbf{X}\|_{F}}\right)$
(3)

where $\textbf{X}\in\mathbb{C}^{n\times m}$ is the desired data matrix,
$\textbf{X}^{\sharp}\in\mathbb{C}^{n\times m}$ is the reconstructed
approximation.

###  III-A Relationship between the SG and quality of reconstruction

We begin with examples that numerically demonstrate the relationship between
the SG and the quality of LRMR output. We perform experiments on 5D data
tensors (corresponding to 3D seismic surveys respectively). We consider the BG
3D model with 68 sources and 401 receivers along each axis to generate 5D data
tensor with dimensions: time, sourcex, sourcey, receiverx, and receivery. We
work in the Fourier domain by applying a fast Fourier transform (FFT) along
the time axis and reconstruct in a frequency by frequency basis (generating a
sequence of 4D tensors). We use the matricization from [17] (first introduced
in [16]) to reorganize each 4D tensor into a matrix, where the sourcex-
receiverx pairs are coupled along the rows and the sourcey-receivery pairs are
coupled along the columns (with the sources serving as the outer dimensions).
Henceforth, we will refer to this matricization as _src-rec form_. We restrict
ourselves to a single 7.34 Hz frequency slice and the first 5 sources along
each axis, generating a $2005\times 2005$ matrix (see Figure 2a).

We begin by removing a specified number of receivers (according to a sampling
percentage) in a periodic manner (e.g., keep only first receiver of every
three for approximately $66\%$ missing receivers). We then solve the nuclear
norm minimization problem (2) with $\eta=0$, applying the corresponding
sampling mask and compute its SG ratio
$\frac{\sigma_{2}(\textbf{M})}{\sigma_{1}(\textbf{M})}$. Subsequently, we
choose a specified percentage $p\in[0,1]$ of the observed receivers and
relocate them among the unobserved receiver locations in a uniform random
manner. This generates sampling masks with the specified sampling percentage
and varying levels of randomly placed receivers. We now solve the nuclear norm
minimization problem with the modified sampling mask and record the respective
SG ratio. We repeat this procedure for several percentage $p$ values in an
increasing manner (specifically, $p\in\\{0,.02,.04.,\cdots,1\\}$).

(a) a

(b) b

(c) c

Figure 2: Examples of recovered data matrices via sampling masks with small
and large SG. (a) True data (b) reconstructed data with sampling mask
exhibiting SG ratio $=.9828$ and recovery SNR $=3.5$ dB ($75\%$ missing
receivers), and (c) reconstructed data with sampling mask exhibiting SG ratio
$=.1796$ with SNR $=20.7$ dB ($75\%$ missing receivers).

In general, we observe a correlation between the number of randomly placed
receivers and the SG of the resulting sampling mask, where masks with larger
degrees of randomness provide larger SG values (i.e., smaller SG ratios
$\frac{\sigma_{2}}{\sigma_{1}}$). This whole procedure is repeated 100 times
independently and we average the SG ratios and SNR values for each
$p\in\\{0,.02,.04.,\cdots,1\\}$. To obtain Figure 1 we sort the resulting SG
ratios in increasing order and plot them against the corresponding SNR of
reconstruction. Furthermore, Figures 2b and 2c show two recovered data sets
from sampling masks with small and large SG, respectively. These figures imply
a clear correlation between the SG and the SNR of reconstruction.

Notice that the SNR values in Figure 1 taper off rather quickly towards the
optimal SNR value in each case. This observation implies that a practitioner
need not place all receivers randomly and might achieve a similar quality of
reconstruction with a hybrid layout (e.g., half periodic and half random
receiver placement). Such acquisition designs provide achievable surveys with
near-optimal SG ratios (i.e., close to those of a Ramanujan graph [13]), where
the resulting survey minimizes the amount of receivers not placed in a
convenient and patterned fashion.

###  III-B Sensitivity of the SG to large gaps of missing data

The previous section implicitly demonstrated the utility of randomness for
LRMR appropriate sampling via the SG. However, such random samples may exhibit
large gaps of missing data that incur degrading artifacts on many
interpolation techniques [38]. In this section, we explore the sensitivity of
the SG to such gaps. For this purpose we use jittered subsampling [12], i.e.,
we partition the desired dense receiver grid into intervals with equal number
of receivers and randomly choose a single receiver location to sample from
each interval (where the number of receivers per interval depends on the
sampling ratio). These sampling techniques are common in computer graphics to
attenuate aliasing artifacts while controlling the gap size of unrendered data
[8].

We design an abstract 3D survey (corresponding to a 5D data tensor) with 32
sources and receivers in each axis and restrict ourselves to a single
frequency slice, working in Fourier domain and src-rec form as in the previous
section (generating a $1024\times 1024$ matrix). To emphasize that our
methodology does not require numerically expensive data simulation, we do not
use an explicit dataset here and instead provide our argument via the binary
sampling masks only.

We remove $80\%$ of the receivers via several jitter subsampling methods that
generate varying average gap sizes of missing receiver data (here each
interval contains 5 receivers). We do so by varying the probability
distribution by which we choose the receiver location for every other
interval, while the remaining intervals retain the uniform distribution (i.e.,
each receiver location has equal probability of being chosen). The selection
distribution is modified by selecting a jitter parameter $:=\rho\in[0,1]$ and
only choosing a sample among the first $\lceil\rho\cdot 5\rceil$ receiver
locations in the interval. This provides receiver locations with increasing
average gap sizes (as $\rho$ decreases). Specifically we consider
$\rho\in\\{.2,.4,.6,.8,1\\}$. By only modifying every other interval’s
selection distribution we still assure a rich sampling scheme (i.e., each row
and column will be sampled with high probability in the chosen matricization).

We compute the SG ratio of each sampling mask generated for each $\rho$ value
specified above. We repeat the experiments 200 times and plot the average SG
ratio against the corresponding jitter parameter $\rho$ in Figure 3. The
experiments are repeated for $90\%$ and $95\%$ missing receivers. Figure 3
exhibits a clear trend between the jitter parameter $\rho$ and the SG, where
$\rho\approx 1$ provides the largest SG (smallest SG ratios
$\frac{\sigma_{2}}{\sigma_{1}}$). The case $\rho=1$ corresponds to optimally-
jittered subsampling [12], with smallest average gap size between receiver
locations. Arguably, this type of jitter is best for general sampling and the
SG can successfully make this distinction. The main conclusion of this section
is that the SG can be used to detect a well distributed set of samples (which
we generate in a random fashion for convenience) while remaining aware of
large gaps of missing data. These properties are crucial for universal
sampling schemes and the SG is an efficient tool to quantify these in a given
set of observed matrix entries.

(a)  Figure 3: Plot of average SG ratio of sampling mask vs jitter parameter
$\rho$, for various sampling percentages. The results are taken from the
average of 200 independent experiments. Notice that as the average sampling
gap increases ($\rho\rightarrow 0$) we obtain less favorable SG values (and
diminished quality of reconstruction as expected). The SG is therefore
sensitive to large gaps of missing data and we expect the highest quality of
reconstruction via optimally-jittered subsampling ($\rho=1$).

###  III-C Coil Sampling Experiments

We consider a 3D marine streamer coil survey, which is designed to have a
random distribution of sources and receivers [26]. For imaging and post-
processing, such data sets must be accurately interpolated onto a dense
regular grid. We focus on this reconstruction step for a dual-coil streamer
data set provided by Schlumberger, deployed on a survey area of $10\times 10$
km (illustrated via src-rec form in Figure 4). We consider the source-receiver
layout indicated by the streamer data, with the goal of choosing an output
interpolation grid that is as dense as possible while appropriate for LRMR.
The experiments of this section do not utilize the actual pressure wave
recordings of the survey (this data remains the property of Schlumberger).
Instead, we verify the quality of reconstruction by simulating seismic data
using the SEG/EAGE Overthrust model that matches the $10\times 10$ km area
from which the coil sampling geometry was extracted.

We adopt the same work flow from the previous sections, reconstructing
sequentially on 4D mono-chromatic tensors matricized via src-rec form and
present our results for a single 10 Hz frequency slice. Given a desired
spacing of sources and receivers, we create a corresponding sampling mask
$\textbf{M}\in\\{0,1\\}^{n\times m}$, with one’s in the entries where streamer
data was acquired and zero’s assigned to all other entries. Note that, in the
field acquisition, sources and receivers are placed off-the-grid (i.e., not
necessarily in our desired output grid). Therefore, to construct the sampling
mask we _bin_ the data onto the desired reconstruction grid, i.e., we assign
each off-the-grid sample to its nearest grid point. We adopt this binning
procedure for simplicity in acquisition design, but note that LRMR
methodologies exist that incorporate the off-the-grid samples for improved
quality of reconstruction [21].

We examined many output grids by varying the source-receiver density between
50-200 m and showcase three examples in Table I along with the respective SG
ratios and SNR of reconstruction (computed for the synthetic 10 Hz Overthrust
frequency slice). The SG successfully distinguishes which layout will produce
the highest quality output in terms of the SNR. Notice that a denser grid is
not always optimal since the example in Table I with densest 50 m source-
receiver spacing exhibits the least desirable SG ratio (and reconstruction
SNR). Intuitively, this 50 m grid is relatively too dense for the given coil
samples and therefore decreases the sampling percentage with large gaps of
missing data that are conveniently detected via the SG.

(a)  Figure 4: Sampling mask generated from a coil-based marine seismic survey
deployed by Schlumberger. The data corresponds to single-coil survey,
extracted from a dual-coil survey. The array represents the entire desired
densely sampled seismic data reshaped as a matrix in src-rec form (see Section
III-A), where each entry corresponds to a source-receiver pair. Blue entries
represent the sampled data, i.e., the receiverx-receivery pairs that recorded
the pressure wave produced by a sourcex-sourcey pair.  TABLE I: SG & SNR AS A
FUNCTION OF GRID SPACING

Receiver spacing (m) | Source spacing (m) | SG ratio ($\frac{\sigma_{2}}{\sigma_{1}}$) | SNR
---|---|---|---
100 | 200 | 0.77 | 10 dB
50 | 100 | 0.4 | 16 dB
50 | 50 | 0.88 | 8 dB

Reconstruction experiments for marine seismic data, where the SG ratios
correspond to sampling masks in src-rec form extracted from a dual-coil survey
deployed by Schlumberger. The SNR corresponds to LRMR experiments using a
simulated 10 hz frequency slice via the SEG/EAGE Overthrust model that matches
the survey area.

###  III-D SG for sparsity-based reconstruction

Although our focus thus far has been on the matrix completion approach, graph
spectrum-based techniques are flexible and apply to other signal processing
and imaging frameworks (including compressive sensing [22]). To validate the
SG on sparsity-based seismic data reconstruction, we designed a 2D survey
layout for an ocean-bottom node (OBN) acquisition and perform reconstruction
using matching pursuit Fourier interpolation [37]. We use synthetic data
simulated on the geologically complex SEAM model. The full dataset contains
510 sources and receivers sampled at 10 m. To simulate a realistic subsampling
scenario, we sample sources at 50 and 100 m across inline and crossline
directions, respectively. This subsampling will cause aliasing to appear from
15 and 7 Hz onwards across the frequency band in the inline and crossline
directions. We only perform reconstruction for a common-receiver gather. We
interpolate a periodic and a jittered random acquisition design [12],
comparing the resulting SG and SNR of reconstruction.

The periodic and jittered subsampling masks have respective SG ratios
$\frac{\sigma_{2}(\textbf{M})}{\sigma_{1}(\textbf{M})}$ of 1 and 0.35. For
periodic sampling, we subsample data at every $5^{th}$ and $10^{th}$ grid
point along inline and crossline directions. Figure 5a shows the cross-line
section from the SEAM dataset, whereas Figures 5b and 5c show the subsampled
data using periodic and jittered sampling respectively. Figures 6a and 6b show
the interpolation results for periodic and jittered subsampling with an SNR
values of 20 and 24 dB, respectively. To stabilize the reconstruction, priors
were introduced in the data reconstruction framework. The amplitude spectrum
derived from the alias-free frequency band of small spatio-temporal windows is
incorporated as weights to distinguish between the aliasing artifacts and the
true events at the higher frequencies [36, 30]. Although the priors
significantly improve the results for the periodic sampled data, we still
loose coherent energy of complex seismic events. In comparison, the
reconstruction via jittered subsampling with priors is able to preserve the
continuity of reflection events and complex diffraction patterns. This
observation is predicted by the SG and further supported by the residuals in
Figures 6c and 6d.

(a)  Figure 5: (a) Fully sampled crossline section extracted from the SEAM
model. Subsampled crosssection using (b) periodic, and (c) jittered
subsampling designs. The SG ratios
($\frac{\sigma_{2}(\textbf{M})}{\sigma_{1}(\textbf{M})}$) for (b) and (c) are
1.0 and 0.35, respectivelty. (a)  Figure 6: Crossline section from the SEAM
model after seismic data reconstruction using (a) periodic, and (b) jittered
subsampling, and (c, d) the corresponding residual sections. The
reconstruction SNRs for (b) and (c) are 20 dB and 24 dB respectively. In
contrast to periodic sampling, we are able to preserve the reflections and
diffraction energy well using the jittered subsampling, which is evident in
the residual section and the SG.

##  IV Conclusions and Discussion

This paper provides a computationally efficient tool that extends the utility
of matrix completion and compressive sensing for seismic data acquisition. By
considering the work in [3, 11, 32], we propose the spectral gap as a means to
select amongst source-receiver layouts in terms of their suitability for
reconstruction via low rank matrix recovery and compressive sensing. Our
numerical experiments demonstrate the flexibility of this tool and its clear
correlation with quality of reconstruction.

The main advantage of the spectral gap is its relatively low computational
complexity, were a practitioner need only compute the first two singular
values of a binary matrix. However, the spectral gap is oblivious to the
underlying data and does not consider many important aspects of seismic survey
design (e.g., azimuth coverage and fold). Because of this, we warn the reader
that the spectral gap does not provide a comprehensive survey design by
itself. Instead, as demonstrated in this work, the spectral gap can be
incorporated into standard survey design workflows. When utilized prior to in
field surveys and large-scale optimization, the spectral gap is a flexible
tool that severely reduces the amount of simulation-based evaluations required
for acquisition design.

Several directions for future work remain. The main results in [3, 11, 32]
require that $\Omega$ be generated from a regular graph (this is the condition
on the singular vectors of the sampling mask M in Theorem II.1), which may be
too restrictive. While our experiments have illustrated that the spectral gap
is informative even when this condition may not hold, a general result that
relaxes this requisite would be most instructive for cost-effective and low-
environmental impact acquisition design.

Throughout, we implicitly assume that the samples are on-the-grid since our
sampling masks are generated by choosing sources-receivers from a dense
equispaced grid. However, in practice this is rarely the case since samples
often deviate from the desired grid points and in general lie on a continuous
domain (as in the dual-coil survey of Section III-C). Therefore, coil sampling
and other randomized acquisition designs would benefit from an analogous tool
that quantifies these off-the-grid source-receiver layouts.

This work has shown that the spectral gap can be a valuable tool when choosing
among a finite set of acquisition designs. Although this alleviates survey
design, practitioners must exhaustively consider many abstract sampling
schemes to be compared via the spectral gap. An improved tool could utilize
the concepts of this work to choose among an infinite set of layouts. This
could be done, for example, by setting up a tractable optimization program
that finds the sparsest mask with largest spectral gap among all masks that
satisfy constraints specified by a desired sampling geometry.

Beyond seismology, one interesting application of the spectral gap could be
for designing new 2D medical ultrasound probes. The total number of
transducers in a probe is a function of the sampling rate based on Nyquist
criteria and the required size of the probe. Technical challenges related to
physical wired connections of a large number of transducers makes it necessary
to reduce this number. This reduction is non trivial, since the quantity of
transducers should not affect the quality of the images generated from probe
measurements. Random distribution of transducers is beneficial for improving
the ultrasound images [9]. As the modelling of a large number of transducers
is time consuming, the spectral gap could speed up the probe design process.

##  V Acknowledgments

The authors would like to thank Schlumberger WesternGeco for their support and
for providing of the dual-coil streamer data.

## References

  * [1] Allouche, N., Øzdemir, K., Özbek, A., & Hopperstad, J. F. [2020]. Multimeasurement aliased-noise attenuation for sparse land seismic acquisition using compressed sensing. _Geophysics_ , 85(2), V183-V200.
  * [2] Aleksandr Aravkin, Rajiv Kumar, Hassan Mansour, Ben Recht, and Felix J. Herrmann [2014]. Fast Methods for Denoising Matrix Completion Formulations, with Applications to Robust Seismic Data Interpolation. _SIAM J. Sci. Comput._ , 36(5), S237-S266.
  * [3] Srinadh Bhojanapalli and Prateek Jain [2014], Universal Matrix Completion. _Proceedings of the 31st International Conference on Machine Learning_ , 32(2), 1881-1889.
  * [4] Emmanuel Candés and Benjamin Recht [2009], Exact Matrix Completion via Convex Optimization. _Foundations of Computational Mathematics_ , 9(6), 717-772.
  * [5] Emmanuel Candés and Terence Tao [2010], The Power of Convex Relaxation: Near Optimal Matrix Completion. _IEEE Transactions of Information Theory_ , 56(5), 2053-2080.
  * [6] Emmanuel Candés and Yaniv Plan [2009], Matrix Completion with Noise. _Proceedings of the IEEE_ , 98(6), 925-936.
  * [7] Y. Chen, S. Bhojanapalli, S. Sanghavi, and R. Ward. [2015], Completing any Low-rank Matrix, Provably. _Journal of Machine Learning Research_ , 16, 2999-3034.
  * [8] Christensen, A. Kensler and C. Kilpatrick [2018], Progressive Multi-Jittered Sample Sequences. _Eurographics Symposium on Rendering_ , 37(4).
  * [9] Cobbold, Richard S.C. [2007], Foundation of Biomedical Ultrasound. _Oxford University Press_.
  * [10] D. L. Donoho [2006], Compressed sensing. _IEEE Transactions on Information Theory_ , 52(4), 1289-1306.
  * [11] E. Heiman, G. Schechtman and A. Schraibman [2014], Determinitic Algorithms for Matrix Completion. _Random Structures and Algorithms_ , 45(2), 306–317.
  * [12] Gilles Hennenfent and Felix Herrmann [2008], Simply Denoise: Wavefield Reconstruction via Jittered Undersampling. _Geophysics_ , 73(3), 1MJ-Z46.
  * [13] Shlomo Hoory, Nathan Linial and Avi Wigderson [2006], Expander Graphs and Their Applications. _Bulletin of the American Mathematical Society_ , 43(4), 439-561.
  * [14] Mike Howard [2007], Marine seismic surveys with enhanced azimuth coverage: Lessons in survey design and acquisition. _The Leading Edge_ , 26, 480-493.
  * [15] P. Jain, P. Netrapalli and S. Sanghavi [2013], Low-rank Matrix Completion Using Alternating Minimization. _Geophysics_ , 78, 1SO-Z134.
  * [16] Kinneging, N. A., Budejicky, V., Wapenaar, C. P. A. and Berkhout, A. J. [1989], Efficient 2D and 3D Shot Record Redatuming. _Geophysical Prospecting_ , 37(5), 493-530.
  * [17] Rajiv Kumar, Curt Da Silva, Okan Akalin, Aleksandr Y. Aravkin, Hassan Mansour, Benjamin Recht and Felix Herrmann [2015], Efficient Matrix Completion for Seismic Data Reconstruction. _Geophysics_ , 80(5), 1SO-Z72.
  * [18] Rajiv Kumar, Oscar López, Damek Davis, Aleksandr Y. Aravkin and Felix Herrmann [2017], Beating Level-Set Methods for 5-D Seismic Data Interpolation: A Primal-Dual Alternating Approach. _IEEE Transactions of Computational Imaging_ , 3(2), 264-274.
  * [19] F. Lazebnik and V. A. Ustimenko [1995], Explicit construction of graphs with an arbitrary large girth and of large size. _Discrete Appl. Math._ , 60(3), 275–284.
  * [20] F. Lazebnik and A. J. Woldar [2001], General properties of some families of graphs defined by systems of equations. _J. Graph Theory_ , 38, 65–86.
  * [21] Oscar López, Rajiv Kumar, Özgür Yılmaz and Felix Herrmann [2016], Off-the-grid low-rank matrix recovery and seismic data reconstruction. _IEEE Journal of Selected Topics in Signal Processing_ , 10(4), 658-671.
  * [22] W. Lu, W. Li, W. Zhang, and S. Xia [2019], Expander recovery performance of bipartite graphs with girth greater than 4. _IEEE Trans. Signal Inf. Process. Netw._ , 5(3), 418–427.
  * [23] Jianwei Ma [2013], Three-dimensional Irregular Seismic Data Reconstruction via Low-rank Matrix Completion. _Geophysics_ , 78, 1SO-Z134.
  * [24] Hansruedi Maurer, Andrew Curtis, and David E. Boerner [2010], Recent advances in optimized geophysical survey design. _GEOPHYSICS_ , 75, A177-A194.
  * [25] Julien Meunier [2011]. Seismic Acquisition from Yesterday to Tomorrow. _Society of Exploration Geophysicists_.
  * [26] Nick Moldoveanu [2010]. Random Sampling: A New Strategy for Marine Acquisition. _SEG Technical Program Expanded Abstracts_ , 51-55.
  * [27] Nicolae Moldoveanu, Nikolas Nesladek and Denes Vigh [2020]. Wide-azimuth towed-streamer and large-scale OBN acquisition: A combined solution. _SEG Technical Program Expanded Abstracts_ , 16-20.
  * [28] Mosher, C. C., Li, C., Janiszewski, F. D., Williams, L. S., Carey, T. C., & Ji, Y. [2017]. Operational deployment of compressive sensing systems for seismic data acquisition. _The Leading Edge_ , 36(8), 661-669.
  * [29] V. Oropeza and M. Sacchi [2011], Simultaneous Seismic Data Denoising and Reconstruction via Multichannel Singular Spectrum Analysis. _Geophysics_ , 76, V25-V32.
  * [30] Özbek, A., M. Vassallo, A. K. Özdemir, D. Molteni, and Y. K. Alp [2010], Anti-alias optimal interpolation with priors: 80th Annual International Meeting, SEG, Expanded Abstracts, 3401–3405.
  * [31] C. Peirce and A. J. Day [2002], Ocean-bottom seismograph tomographic experiments—a consideration of acquisition geometries vs. resources. _Geophysical Journal International_ , 151(2), 543-565.
  * [32] Shantanu Prasad Burnwal and Mathukumalli Vidyasagar [2020], Deterministic Completion of Rectangular Matrices Using Asymmetric Ramanujan Graphs: Exact and Stable Recovery. _IEEE Transactions of Signal Processing_ , 68, 3834-3848.
  * [33] Benjamin Recht [2011], A Simpler Approach to Matrix Completion. _Journal of Machine Learning Research_ , 12, 3413-3430.
  * [34] Benjamin Recht and Christopher Ré [2013], Parallel Stochastic Gradient Algorithms for large-scale Matrix Completion. _Mathematical Programming Computation_ , 5, 201-226.
  * [35] Benjamin Recht, Maryam Fazel and Pablo A. Parrilo [2010], Guaranteed Minimum-Rank Solutions of Linear Matrix Equations via Nuclear Norm Minimization. _SIAM Rev._ , 52(3), 471-501.
  * [36] Schonewille, M., A. Klaedtke, A. Vigner, J. Brittan, and T. Martin [2009], Seismic data regularization with the anti-alias antileakage Fourier transform. _First Break_ , 27(9), 85–92.
  * [37] Schonewille, M., Yan, Z., Bayly, M., & Bisley, R. [2013]. Matching pursuit Fourier interpolation using priors derived from a second data set. _Society of Exploration Geophysicists Technical Program Expanded Abstracts_ , 3651-3655.
  * [38] D.O. Trad, J. Deere and S. Cheadle [2005], Challenges for land data interpolation: Presented at the CSEG National Convention.
  * [39] E. J. van Veldhuizen, G Blacquière, and A. J. Berkhout [2008], Acquisition geometry analysis in complex 3D media. _GEOPHYSICS_ , 73, Q43-Q58.
  * [40] Gijs J.O. Vermeer [2012], 3D Seismic Survey Design. _Society of Exploration Geophysicists, Okla., USA_ , Second Edition.
  * [41] S. Wu, D. J. Verschuur and G. Blacquière [2022], Automated Seismic Acquisition Geometry Design for Optimized Illumination at the Target: A Linearized Approach. _IEEE Transactions on Geoscience and Remote Sensing_ , 60, 1-13.
  * [42] Vinje, Vettel and Lie, Jan Erik and Danielsen, Vidar and Dhelle, Per Eivind and Silliqi, Risto and Nilsen, Carl-Inge and Hicks, Erik and Camerer, Anne [2017], Shooting over the seismic spread. _First Break_ , 35(6), 1365-2397.
  * [43] Takamitsu Watanabe and Naoki Masuda [2010], Enhancing the Spectral Gap of Networks by Node Removal. _Phys. Rev. E_ , 82(4).
  * [44] Yi Yang, Jianwei Ma and Stanley Osher [2013], Seismic data reconstruction via matrix completion. _Inverse Problems and Imaging_ , 7(4), 1379-1392.
  * [45] H. Yun, H. Yu, C. Hsieh, S. Vishwanathan, & I. Dhillon [2014], NOMAD: Non-locking, stochastic Multi-machine algorithm for Asynchronous and Decentralized matrix completion. _Proceedings of the VLDB Endowment_.