# Transmission in graphene through tilted barrier in laser field

Rachid El Aitouni  Laboratory of Theoretical Physics, Faculty of Sciences,
Chouaïb Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Miloud
Mekkaoui  Laboratory of Theoretical Physics, Faculty of Sciences, Chouaïb
Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Ahmed Jellal
<EMAIL_ADDRESS>Laboratory of Theoretical Physics, Faculty of Sciences,
Chouaïb Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Canadian
Quantum Research Center, 204-3002 32 Ave Vernon, BC V1T 2L7, Canada

###### Abstract

We study the transmission of Dirac fermions in graphene through a tilted
barrier potential in the presence of a laser field. By using Floquet theory,
we solve the Dirac equation and then obtain the energy spectrum. The boundary
conditions together with the transfer matrix method allow us to determine the
transmission probabilities corresponding to all energy bands. For numerical
purposes, we only study the central band $l=0$ and the two first side bands
$l=\pm 1$ under various conditions. It has been discovered that the laser
field and potential influence the amplitude and shape of the transmission
probabilities. In fact, different oscillations are observed, the Klein effect
is still present, and the transmission of all bands is almost zero inside the
barrier.

###### pacs:

78.67.Wj, 05.40.-a, 05.60.-k, 72.80.Vp
Keywords: Graphene, laser field, barrier potentiel, transmission probability.

##  I Introduction

Graphene is one of the few two-dimensional materials that is made up of carbon
atoms that form a hexagonal structure in the shape of a honeycomb [1]. The
thickness of which does not exceed the thickness of an atom. It has incredible
mechanical, electronic, optical, thermal, and chemical properties [2], as well
as extremely high mobility. The speed of its electrons is 300 times less than
the speed of light [3]. Because of its speed, the researchers consider its
electrons to be massless Dirac fermions. It is one of the most powerful
conductors of electricity and heat, which has started to be used in technology
industries [4, 5, 6]. It replaces ordinary semiconductors such as sulisum to
manufacture electronic components, solar panels, and touch displays. The
Hamiltonian describing the graphene properties is a relativistic Dirac
Hamiltonian, resulting within the framework of the tight-binding model [3].
The eigenvalue equation resolution reveals that the Brillouin zone is
delimited by six Dirac points, which are represented by two non-equivalent
points $K$and $K^{\prime}$, each corresponding to two atoms of the pattern in
the direct lattice. In addition, the energy is a linear dispersion relation in
the vicinity of the Dirac points [7], which is similar to a cone, i.e., the
valence and the conduction band are tangent at the Dirac points. As a result,
the electrons inside graphene can easily jump from the valence band to the
conduction band. This is the problem that delays its use. The zero gap between
these two bands prompted many researchers to look for ways to create a gap
between them. Among the methods proposed are: confining the charge to the
surface in systems composed of multi-layered graphene [8], deformation of the
graphene sheet [3, 9, 10, 11, 12], and application of different fields like
electric, magnetic, or laser fields [13, 14, 15, 16, 17].

Laser technology has become a powerful tool for the advancement of
investigations into graphene. Indeed, several studies have been conducted on
the effect of a linear or circular laser field on the mobility of Dirac
fermions in graphene [18, 19, 20, 18], demonstrating that the laser can create
a band gap in the energy spectrum. On the other hand, the tunneling effect has
been studied for Dirac fermions in graphene subjected to a linear vector
potential [13, 21]. Here, the infinite mass boundary condition reduced our 2D
Dirac equation to a massively effective 1D Dirac equation with an effective
mass equal to the quantized transverse momentum. Note that the transmission
probabilities of Dirac fermions through different barriers are studied. In
particular, rectangular barriers [13] and inclined barrier [21] all show that
the transmissions are not null even if the incident energies are lower than
the barrier height, which is called the Klein’s tunnel effect [22, 23].
Barrier potentials in conjunction with laser fields have also been
investigated [24, 15, 16] and potential oscillations over time [25]. As a
result, it was discovered that side bands $l\hbar\omega\ (l=0,\pm 1,\cdots)$
were added to the energy spectrum, giving rise to an infinite number of
transmission probabilities.

We investigate how Dirac fermions in graphene can pass through a slanted
barrier potential in the presence of a laser field. The Dirac equation is
solved using Floquet theory, and the energy spectrum is then obtained. The
transfer matrix approach and boundary conditions enable us to ascertain the
transmission probabilities corresponding to all energy bands. Under different
conditions, only the central band $l=0$ and the two first side bands $l=\pm 1$
are examined numerically. It has been found that the amplitude and shape of
the transmission probabilities are affected by the laser field and potential.
In reality, a variety of oscillations are seen inside the barrier, the Klien
effect is still evident, and the transmission of all bands is practically
nonexistent.

This paper is organized as follows. In Sec. II, we set the mathematical
formalism based on the Hamiltonian describing the present system. By solving
the eigenvalue equation, we will be able to explicitly determine the
eigenspinors in all regions composing the system under consideration. These
solutions will allow us, in Sec. III, to use the boundary conditions at
interfaces, and by employing the transfer matrix approach, we will calculate
all transmission channels resulting from the oscillating barrier, which causes
the energy subbands $l\hbar\omega\ (l=0,\pm 1,\ldots)$. In Sec. IV, we will
numerically analyze our results by plotting the three first transmission
probabilities as a function of the physical parameters. Finally, we conclude
our work.

##  II Theoretical mode

We consider a graphene sheet subjected to a tilted potential $V(x)$ and
irradiated by a linear polarization monochromatic laser $A(t)$ over a finite
region of length $d$, while the other two regions are pure graphene as
depicted in Fig. 1.

Figure 1:  Shematic diagram of a graphene sheet electrostatic inclined barrier
irradiated by a monochromatic laser field over a finite region of length $d$.

The Hamiltonian describing the system is given by

$H=v_{F}\vec{\sigma}\cdot\left[\vec{p}-e\vec{A}(t)\right]+V(x)\mathbb{I}_{2}$
(1)

where $v_{F}$ is the Fermi velocity $\approx 10^{6}m/s$, $\sigma_{x,y}$ are
the Pauli matrices, $\vec{p}=-i\hbar(\partial x,\partial y)$ is two
dimensional momentum, $e$ being the electronic charge and $\mathbb{I}_{2}$ is
$2\times 2$ identity matrix. The scalar potential $V(x)$ varies along the
$x$-axis and is represented by

$\displaystyle V(x)=\begin{aligned} &\begin{cases}-\beta x+V_{0},&0<x<d\\\
0,&\text{otherwise}\end{cases}\end{aligned}$ (2)

which can be produced by two identical plates located at $x=0$ and $x=d$,
where $\beta=\frac{V_{0}-V_{1}}{d}$ and $V_{0}>V_{1}$. In the dipole
approximation [26], the vector potential $\vec{A}(t)$ of the laser field is
given by

$\displaystyle\vec{A}(\vec{r},t)=(0,A_{0}\cos(\omega t),0,t)$ (3)

where $A_{0}=\frac{F}{\omega}$ is the laser field amplitude. An electric field
$\vec{E}$ with frequency $\omega$ and amplitude $F$ can generate such a
potential. Indeed, we can write

$\vec{E}=-\frac{\partial}{\partial t}\vec{A}(t)=\vec{F}\sin(\omega t).$ (4)

The present system is composed of three regions labeled $j=1,2,3$. Then we
have to determine the solutions of the energy spectrum in each region using
the eigenvalue equation

$\left[v_{F}\sigma_{x}p_{x}+v_{F}\sigma_{y}(p_{y}-eA(t))+V_{j}(x)\mathbb{I}_{2}\right]\Psi_{j}(x,y,t)=i\hbar\frac{\partial}{\partial
t}\Psi_{j}(x,y,t).$ (5)

Since the system has a finite width $R$, then one can use the boundary
conditions at $y=0,R$

$\displaystyle\Psi_{1}(x,0,t)=-\Psi_{2}(x,0,t)$ (6)
$\displaystyle\Psi_{1}(x,R,t)=-\Psi_{2}(x,R,t)$ (7)

to quantize the transverse wave vector

$\displaystyle k_{y}=\frac{\pi}{R}\left(n+\frac{1}{2}\right),\quad n=0,\pm
1,\cdots.$ (8)

The laser field is considered to oscillate periodically in time, which has an
effect on the behavior of the wave function and the energy band. Consequently,
the Floquet approximation [27] should be taken into account to determine the
solutions of the eigenvalue equation. The eigenspinors can also be separated
in coordinates because the Hamiltonian is invariant along the $y$-direction.
Combining all of these, we can write the eigenspinors solution of (5) as

$\Psi_{j}(x,y,t)=(\psi_{j1}(x),\psi_{j2}(x))^{T}e^{ik_{y}y}\phi_{0}(t)e^{-iEt}$
(9)

where $E$ denotes the Floquet energy and $\phi_{0}(t)$ varies over time,
$\phi_{0}(t+\tau)=\phi_{0}(t)$, $\tau$ denotes the laser field period and $t$
denotes transpose. As a result, (5) provides us with

$\displaystyle\left(-i\frac{\partial}{\partial
x}-i\left(k_{y}-\frac{F}{\omega}\cos(\omega
t)\right)\right)\psi_{j2}(x)\phi_{0}(t)=i\psi_{j1}(x)\frac{\partial\phi_{0}(t)}{\partial
t}+(E-V(x))\psi_{j1}(x)\phi_{0}(t)$ (10)
$\displaystyle\left(-i\frac{\partial}{\partial
x}+i\left(k_{y}-\frac{F}{\omega}\cos(\omega
t)\right)\right)\psi_{j1}(x)\phi_{0}(t)=i\psi_{j2}(x)\frac{\partial\phi_{0}(t)}{\partial
t}+(E-V(x))\psi_{j2}(x)\phi_{0}(t)$ (11)

where dimensionless quantities are introduced $A\to A{F}t_{0}$,
$(x,y)\to(x,y)d_{0}$, $k_{y}\to k_{y}/d_{0}$, $V(x)\to V(x)E_{0}$, $E\to
EE_{0}$, $t\to tt_{0}$ and $\omega\to\omega\omega_{0}$, with the length
$d_{0}=\sqrt{\frac{\hbar v_{F}}{eF}}$, energy $E_{0}=\frac{\hbar
v_{F}}{d_{0}}$ and time $t_{0}=\frac{d_{0}}{v_{F}}$ scales.

It is worth noting that there are only two differential equations and three
unknown wave functions, $\psi_{j1}(x)$, $\psi_{j2}(x)$, and $\phi_{0}(t)$. To
overcome this situation, we apply the iterative method as a first
approximation by supposing that $\psi_{21}$ and $\psi_{22}$ satisfy the Dirac
equation in region 2 without laser irradiation. As a result, (10-11) becomes

$\displaystyle\frac{F}{w}\cos(\omega
t)\psi_{22}(x)\phi_{0}(t)=\psi_{21}(x)\frac{\partial}{\partial t}\phi_{0}(t)$
(12) $\displaystyle-\frac{F}{w}\cos(\omega
t)\psi_{21}(x)\phi_{0}(t)=\psi_{22}(x)\frac{\partial}{\partial t}\phi_{0}(t).$
(13)

From these we end up with the second order differential equation

$\frac{\partial^{2}\phi_{0}(t)}{\partial t^{2}}+\omega\tan(\omega
t)\frac{\partial\phi_{0}(t)}{\partial
t}+\frac{F^{2}}{\omega^{2}}\cos^{2}(\omega t)\phi_{0}(t)=0$ (14)

showing the solution [16]

$\displaystyle\phi_{0}(t)\approx e^{-i\frac{F}{\omega^{2}}\sin(\omega t)}.$
(15)

Injecting this into (10) to get

$\displaystyle\left\\{\left(-i\frac{\partial}{\partial
x}-i\left(k_{y}-\frac{\partial}{\partial
t}\right)\right)\psi_{22}(x)-\left(i\frac{\partial}{\partial
t}+E-V(x)\right)\psi_{21}(x)\right\\}e^{-i\frac{F}{\omega^{2}}\sin(\omega
t)}e^{-im\omega t}=0.$ (16)

Generating the Bessel function $J_{m}$ using the Jacobi-Anger expansion

$\displaystyle e^{-i\frac{F}{\omega^{2}}\sin(\omega
t)}=\sum_{m=-\infty}^{+\infty}J_{m}\left(\frac{F}{\omega^{2}}\right)e^{-im\omega
t}$ (17)

and then we show that (16) can be reduced to

$\displaystyle\sum_{m=-\infty}^{+\infty}\left\\{\left(-i\frac{\partial}{\partial
x}-i\left(k_{y}-\frac{\partial}{\partial
t}\right)\right)\psi_{22}(x)-\left(i\frac{\partial}{\partial
t}+E-V(x)\right)\psi_{21}(x)\right\\}J_{m}\left(\frac{F}{\omega^{2}}\right)e^{-im\omega
t}=0.$ (18)

Consequently, we end up with

$\displaystyle\left(V(x)-(E+l\omega)\right)\psi_{21}(x)=-i\left(\frac{\partial}{\partial
x}+(k_{y}-l\omega)\right)\psi_{22}(x).$ (19)

Similarly, we get from (11)

$\displaystyle\left(V(x)-(E+l\omega)\right)\psi_{22}(x)=-i\left(\frac{\partial}{\partial
x}-(k_{y}-l\omega)\right)\psi_{21}(x)$ (20)

Let us make a change to the spinor components to obtain the Schrödinger
equation. This is

$\displaystyle\psi_{21}(x)=\chi_{1}+i\chi_{2}$ (21)
$\displaystyle\psi_{22}(x)=\chi_{1}-i\chi_{2}$ (22)

After substitution in (19-20), we obtain

$\frac{\partial^{2}\chi_{1}}{\partial x^{2}}+\left(-i\frac{\partial
V(x)}{\partial
x}+(k_{y}-l\omega)^{2}+\left(V(x)-(E+l\omega)\right)^{2}\right)\chi_{1}=0$
(23)

which can be written as Weber’s differential equation [28]

$\frac{\partial^{2}\chi_{1}(z_{l})}{\partial{z_{l}}^{2}}+\left(\frac{1}{2}-v_{l}-\frac{z^{2}_{l}}{4}\right)\chi_{1}(z_{l})=0$
(24)

where we have set

$\displaystyle z_{l}=\sqrt{\frac{2}{\beta}}e^{\frac{i\pi}{4}}(\beta x+E_{l})$
(25) $\displaystyle v_{l}=\frac{i}{2\beta}\left(k_{y}-l\omega\right)^{2}$ (26)
$\displaystyle E_{l}=E+l\omega-V_{0}.$ (27)

The solution of (24) is given by

$\displaystyle\chi_{1}(z_{l})=C_{1,l}D_{v_{l}-1}(z_{l})+C_{2,l}D_{-v_{l}}(-z^{*}_{l})$
(28)

where $C_{1,l}$ and $C_{2,l}$ are two constants. From $\chi_{1}$, we derive
the second component

$\displaystyle\chi_{2}(z_{l})=C_{1,l}\frac{\sqrt{2\beta}e^{-i\frac{\pi}{4}}}{k_{y}-l\omega}D_{v_{l}}(z_{l})+\frac{C_{2,l}}{(k_{y}-l\omega)}\left(2\left(\beta
x+E_{l}\right)D_{-v_{l}}(-z^{*}_{l})+\frac{\sqrt{2B}e^{i\frac{\pi}{4}}}{k_{y}-l\omega}D_{1-v_{l}}(-z^{*}_{l})\right).$
(29)

These allow us to express $\psi_{21}(x)$ and $\psi_{22}(x)$ as follows

$\displaystyle\psi_{21}(x)=C_{1,l}\mu^{+}_{l}(x)+C_{2,l}\zeta^{+}_{l}(x)$ (30)
$\displaystyle\psi_{22}(x)=C_{1,l}\mu^{-}(x)+C_{2,l}\zeta^{-}(x)$ (31)

where we shave defined

$\displaystyle\mu_{l}^{\pm}(x)$ $\displaystyle=$ $\displaystyle
D_{v_{l}-1}(z_{l})\pm\frac{\sqrt{2\beta}}{k_{y}-l\omega}e^{\frac{i\pi}{4}}D_{v_{l}}(z_{l})$
(32) $\displaystyle\zeta_{l}^{\pm}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{k_{y}-l\omega}\left(\left(k_{y}-l\omega\mp(2i\beta
x+2iE_{l})\right)D_{-v_{l}}(-z^{*}_{l})\pm\sqrt{2\beta}e^{-i\frac{\pi}{4}}D_{1-v_{l}}(-z^{*}_{l})\right).$
(33)

Combining all to write the eigenspinors in region 2 as

$\Psi_{2}(x,y,t)=e^{ik_{y}y}\sum_{l=-\infty}^{+\infty}\left[C_{1,l}\begin{pmatrix}\mu^{+}_{l}(x)\\\
\ \mu^{-}_{l}(x)\end{pmatrix}+C_{l,2}\begin{pmatrix}\zeta^{+}_{l}(x)\\\
\zeta^{-}_{l}(x)\end{pmatrix}\right]\sum_{m=-\infty}^{+\infty}J_{m}\left(\frac{F}{\omega^{2}}\right)e^{-i(E+l\omega+m\omega)t}.$
(34)

In the regions 1 and 3 have only pristine graphene, the spinors in both
regions written as follows [21]

$\displaystyle\Psi_{1}(x,y,t)=e^{ik_{y}y}\left[\delta_{l,0}\begin{pmatrix}1\\\
\Lambda_{l}\end{pmatrix}e^{ik_{l}x}e^{-iEt}+\sum_{l=-\infty}^{+\infty}r_{l}\begin{pmatrix}1\\\
-\Lambda^{*}_{l}\end{pmatrix}e^{-ik_{l}x}e^{-i(E+l\omega)t}\right]\sum_{m=-\infty}^{+\infty}J_{m}\left(\frac{F}{\omega^{2}}\right)e^{-im\omega
t}$ (35)
$\displaystyle\Psi_{3}(x,y,t)=e^{ik_{y}y}\sum_{m=-\infty}^{+\infty}\left[t_{l}\begin{pmatrix}1\\\
\Lambda_{l}\end{pmatrix}e^{ik_{l}x}+b_{l}\begin{pmatrix}1\\\
-\Lambda^{*}_{l}\end{pmatrix}e^{-ik_{l}x}\right]\delta_{m,l}e^{-i(E+m\omega)t}$
(36)

corresponding to the eigenvalues

$\displaystyle E+l\omega=s_{l}\sqrt{k^{2}_{l}+k^{2}_{y}}$ (37)

with the reflection $r_{l}$ and transmission $t_{l}$ amplitudes,
$\Lambda_{l}=s_{l}e^{i\phi_{l}}$, $\phi_{l}=\tan^{-1}\frac{k_{y}}{k_{l}}$,
$\delta_{m,l}=J_{m-l}(0)$, $s_{l}=\text{sgn}(E+l\omega)$ and ${b_{l}}$ is null
vector.

##  III Transmission probabilities

To determine the transmission and reflection amplitudes, we use the boundary
conditions at the interfaces $x=0$ and $x=d$

$\displaystyle\Psi_{1}(0,y,t)=\Psi_{2}(0,y,t)$ (38)
$\displaystyle\Psi_{2}(d,y,t)=\Psi_{3}(d,y,t)$ (39)

together with the orthogonality of $e^{im\omega t}$ to end up with four
unknown parameters, each one there in an infinity of modes

$\displaystyle\delta_{m,0}+r_{m}=\sum_{l=-\infty}^{+\infty}\left(c_{1,l}\mu^{+}_{l}(0)+c_{2,l}\zeta^{+}_{l}(0)\right)J_{m-l}\left(\frac{F}{\omega^{2}}\right)$
(40)
$\displaystyle\delta_{m,0}\Lambda_{m}-r_{m}\Lambda_{m}^{*}=\sum_{l=-\infty}^{+\infty}\left(c_{1,l}\mu^{-}_{l}(x)(0)+c_{2,l}\zeta^{-}_{l}(0)\right)J_{m-l}\left(\frac{F}{\omega^{2}}\right)$
(41) $\displaystyle
t_{m}e^{ik_{m}d}+b_{m}e^{-ik_{m}d}=\sum_{l=-\infty}^{+\infty}\left(c_{1,l}\mu^{+}_{l}(d)+c_{2,l}\zeta^{+}_{l}(d)\right)J_{m-l}\left(\frac{F}{\omega^{2}}\right)$
(42) $\displaystyle
t_{m}\Lambda_{m}e^{ik_{m}d}-b_{m}\Lambda_{m}^{*}e^{-ik_{m}d}=\sum_{l=-\infty}^{+\infty}\left(c_{1,l}\mu^{-}_{l}(d)+c_{2,l}\zeta^{-}_{l}(d)\right)J_{m-l}\left(\frac{F}{\omega^{2}}\right)$
(43)

These relations can be expressed using the transfer matrix approach to connect
the three regions of the system under consideration. This is

$\begin{pmatrix}\Gamma_{1}\\\
\Gamma^{\prime}_{1}\end{pmatrix}=\mathbb{M}\begin{pmatrix}\Gamma_{2}\\\
\Gamma^{\prime}_{2}\end{pmatrix}$ (44)

with the total transfer matrix

$\displaystyle\mathbb{M}=\begin{pmatrix}M_{1,1}&M_{1,2}\\\
M_{2,1}&M_{2,2}\end{pmatrix}=\mathbb{M}(0,1).\mathbb{M}(1,2)$ (45)

where $\mathbb{M}(j,j+1)$ is the transfer matrix connecting regions $j$ and
$j+1$

$\displaystyle\mathbb{M}(0,1)=\begin{pmatrix}\mathbb{I}&\mathbb{I}\\\
\mathbb{S}^{+}&\mathbb{S}^{-}\\\
\end{pmatrix}^{-1}\begin{pmatrix}\mathbb{F}^{+}_{0}&\mathbb{R}^{+}_{0}\\\
\mathbb{F}^{-}_{0}&\mathbb{R}^{-}_{0}\end{pmatrix}$ (46)
$\displaystyle\mathbb{M}(1,2)=\begin{pmatrix}\mathbb{F}^{+}_{d}&\mathbb{R}^{+}_{d}\\\
\mathbb{F}^{-}_{d}&\mathbb{R}^{-}_{d}\end{pmatrix}^{-1}\begin{pmatrix}\mathbb{I}&\mathbb{I}\\\
\mathbb{S}^{+}&\mathbb{S}^{-}\\\
\end{pmatrix}\begin{pmatrix}\mathbb{Q}^{+}&\mathbb{O}\\\
\mathbb{O}&\mathbb{Q}^{-}\\\ \end{pmatrix}$ (47)

and the quantities

$\displaystyle\mathbb{S}^{\pm}=\pm\delta_{m,l}\Lambda_{l}^{\pm
1},\quad\mathbb{F}^{\pm}_{z}=\mu_{l}^{\pm}(z)J_{m-l}\left(\frac{F}{\omega^{2}}\right)$
(48)
$\displaystyle\mathbb{R}^{\pm}_{z}=\zeta_{l}^{\pm}(z)J_{m-l}\left(\frac{F}{\omega^{2}}\right),\quad\mathbb{Q}^{\pm}=e^{\pm
ik_{l}L}\delta_{m,l}$ (49)

where $\mathbb{O}$ is the zero matrix, $\mathbb{I}$ is the unit matrix and
$z=\\{0,d\\}$. We consider the fermions moving from left to right with the
energy $E$, then from Eq. (44) we get

$\Gamma_{2}=\mathbb{M}^{-1}_{1,1}\Gamma_{1}$ (50)

where $\Gamma_{1}=\delta_{0,l}$ is the Kronecker symbol and
$\Gamma_{2}=t_{l}$. Because $n,l$ varies from $-\infty$ to $+\infty$, the
transfer matrix approach, which we determined above, is of infinite order,
making it difficult to solve. To accomplish this, we replace the infinite
series with a finite number of terms ranging from $-N$ to $N$, where
$N>\frac{F}{\omega^{2}}$, and we limit ourselves only to the first energy side
bands $E\pm m\omega$, yielding the result

$t_{-N+k}=M^{\prime}[k+1,N+1]$ (51)

with $M^{\prime}=M^{-1}_{1,1}$ and $k=0,1,2,\cdots N$.

On the other hand, the transmission probability can be expressed by using the
continuity equation. According to our Hamiltonian, the general expression of
the current density takes the following form:

$\displaystyle J=\Psi_{j}^{*}(x,y,t)\sigma_{x}\Psi_{j}(x,y,t)$ (52)

which makes it possible to determine the current densities in each region $j$.
There are given by

$\displaystyle J_{in,0}=(\Lambda_{0}+\Lambda^{*}_{0})$ (53) $\displaystyle
J_{re,l}=r^{*}_{l}r_{l}(\Lambda_{l}+\Lambda^{*}_{l})$ (54) $\displaystyle
J_{tr,l}=t^{*}_{l}t_{l}(\Lambda_{l}+\Lambda^{*}_{l}).$ (55)

Then by using the relations

$T_{l}=\frac{|J_{tr,l}|}{|J_{in,0}|},\quad
R_{l}=\frac{|J_{re,l}|}{|J_{in,0}|}$ (56)

one can express of the transmission $T_{l}$ and reflection $R_{l}$
probabilities as

$\displaystyle
T_{l}=s^{\prime}\frac{\cos{\theta_{l}}}{\cos{\theta_{0}}}|t_{l}|^{2},\quad
R_{l}=s^{\prime}\frac{\cos{\theta_{l}}}{\cos{\theta_{0}}}|r_{l}|^{2}$ (57)

where the two angles are

$\displaystyle\cos{\theta_{l}}=\frac{k_{l}}{\sqrt{k_{l}^{2}+k_{y}^{2}}},\quad\cos{\theta_{0}}=\frac{k_{0}}{\sqrt{k_{0}^{2}+k_{y}^{2}}}$
(58)

and the sign is $s^{\prime}=s_{l}s_{0}$.

To simplify the numerical calculation, we pronounce only the first sidebands,
i.e. $l=-1,0,1$. The corresponding transmissions are

$\displaystyle t_{-1}=M^{\prime}[1,2]$ (59) $\displaystyle
t_{0}=M^{\prime}[2,2]$ (60) $\displaystyle t_{1}=M^{\prime}[3,2].$ (61)

In the following section, we will plot the transmission probability as a
function of the current system’s physical parameters for the three first
bands: the central band $(l=0)$ and the two first side bands ($l=\pm 1$). We
will analyze them to demonstrate the effect of the barrier that we applied on
the movement of the Dirac fermions.

##  IV Numerical results

Fig. 2 represents the transmission probabilities as a function of the incident
energy $E$. The total transmission $T=T_{0}$ (magenta line) is shown in Fig.
2a, and the other bands $T_{\pm 1}$ are zero in the absence of laser
irradiation $F=0$ of a linear barrier ($V_{0}=20$, $V_{1}=0$). Let us now
examine what happens in the intermediate region if laser irradiation $F=0.9$
is introduced, as shown by the transmission for the central band $T_{0}$ (blue
line) and the two first sidebands $T_{-1}$ (red line) and $T_{1}$ (green
line). We also find that the laser irradiation is very significant in
determining the relative transmission probabilities of two first sidebands.
Transmissions are prohibited for the energy color red, $E=\pm\omega$,
resulting in a change in $k_{y}$ to $k_{y}\pm\omega$, which serves as an
effective mass $m^{*}$ [29, 25]. The total transmission and the transmission
of the central band oscillate around a minimum for the energies $EV_{0}$, then
increase passing by $EV_{0}$ up to a maximum value, and the total transmission
(the sum of the three) tends towards 1, indicating that our system exhibits
classic behavior after the barrier, indicating the presence of the Klein
effect. The same transmissions are plotted in Fig. 2b, but with a two-
potential barrier $(V_{0}=20,V_{1}=15)$. The transmission of the central band
is more dominant in the energy zone $k_{y}<E<V_{1}$, and the total
transmission oscillates in the neighborhood of $1$. The Klein paradox is
always present. We observe a total reflection in the energy zone
$V_{1}<E<V_{0}$, with a very small transmission for the band corresponding to
$l=1$. In the third zone, $E>V_{1}$, total transmission rapidly decreases
towards unit, and $T_{0}$ increases with an apparent oscillation. Finally, we
see that changing the barrier’s inclinison coefficient,
$\beta=\frac{V_{0}-V_{1}}{d}$, affects the shape of the transmissions and the
number of oscillations.

(a)

(b)

Figure 2: (Color online) Transmission probabilities as a function of the
incident energy $E$ with $d=7$, $V_{0}=20$, $k_{y}=1$, $\omega=1.4$, $F=0$ for
$T$ (magenta line), $F=0.9$ for the central transmission band $T_{0}$ (blue
line) and two first sidebands $T_{-1}$ (red line), $T_{1}$ (green line). (a):
$V_{1}=0$ and (b): $V_{1}=15$.

Fig. 3 illustrates the transmission of the central band ($l=0$) as a function
of the incident energy $E$ with $d=7$, $V_{0}=20$, $V_{1}=15$, $k_{y}=1$ for
$\omega=1.4$, $F=\\{0.2,0.9,1.8\\}$ in Fig. 3a and $F=1.8$,
$\omega=\\{1.5,2,2.5\\}$ in Fig. 3b. They show that the transmission
probability varies in an oscillatory way, and guards even allure, but with an
attenuation proportional to $F$ and $\omega$, and oscillations appear after
the barrier. The oscillation becomes important when $F$ increases, as shown in
Fig. 3a, but it also becomes important when $\omega$ decreases, as shown in
Fig. 3b. Finally, the ratio $\frac{F}{\omega^{2}}$ is important because it
determines the shape, amplitude, and oscillations of the transmission
probability. Inside the barrier, we always have a total reflection for all
values of $\omega$ and $F$.

(a)

(b)

Figure 3: (Color online) Transmission probability for the central band $T_{0}$
as a function of the incident energy $E$ with $d=7$, $V_{0}=20$, $V_{1}=15$,
$k_{y}=1$. (a): $\omega=1.4$, $F=0.2$ (green color), $F=0.9$ (red line) and
$F=1.8$ (blue line). (b): $F=1.8$, $\omega=1.5$ (green color), $\omega=2$ (red
line) and $\omega=2.5$ (blue line).

Fig. 4 depicts transmissions as a function of the barrier potential $V_{1}$
under various conditions. Taking $\omega=1.4$ and $F=0.9$ in Fig. 4a, we can
see that transmissions decrease from zero to $E-2k_{y}$, then increase to
$E+2k_{y}$. When $V_{1}>E+2k_{y}$, the transmissions are maximized, and we
have a total transmission. Fig. 4b depicts the transmission $T_{1}$ as a
function of $V_{1}$ with varying $\omega$ and $F$ values. We observe that when
the ratio $\frac{F}{\omega^{2}}$ decreases, the amplitude of $T_{1}$ also
decreases. The green line corresponds to $\frac{F}{\omega}=0.91$ and the blue
line corresponds to $\frac{F}{\omega^{2}}=0.2$. As a result, we note that the
ratio $\frac{F}{\omega^{2}}$ plays an important role in increasing or
decreasing the amplitude of the transmissions.

(a)

(b)

Figure 4: (Color online) Transmission probabilities as a function of potential
$V_{1}$ with $d=7$, $V_{0}=40$, $E=20$, $k_{y}=1$. (a): $w=1.4$, $F=0.9$ for
the total transmission $T$ (magenta line), the central transmission band
$T_{0}$ (blue line) and two first sidebands $T_{-1}$ (red line), $T_{1}$
(green line). (b): Transmission probability $T_{1}$ with
$\\{F=0.4,\omega=1.4\\}$ (blue line), $\\{F=0.9,\omega=1.4\\}$ (red line),
$\\{F=1.8,\omega=1.4\\}$ (green line), $\\{F=1.8,\omega=1.55\\}$ (magenta
line), $\\{F=1.8,\omega=1.7\\}$ (orange line).

Fig. 5 shows the transmission probabilities as a function of the barrier
potential $V_{0}$ with $d=7$, $V_{1}=0$, $E=20$ and $k_{y}=1$. When we use
$\omega=1.4$ and $\frac{F}{\omega^{2}}=0.46$ in Fig. 5a, we see that $T_{0}$
is always more dominant, and $T$ for $V_{1}<E-2k_{y}$ is nearly equal to $1$
(Klein zone). Both $T_{0}$ and $T$ decrease from $E-2k_{y}$ to $E+2k_{y}$,
then increase in an oscillatory way. In Fig. 5b, we only show $T_{-1}$ as a
function of $V_{1}$ with different $\omega$ values. The amplitude and shape of
$T_{-1}$ are clearly dependent on the parameters $\omega$ and $F$. Indeed, the
amplitude is maximum and greater than $0.25$ for $\frac{F}{\omega^{2}}=0.91$
(blue line), then it rapidly decreases to zero after increasing in an
oscillatory manner. The amplitude increases up to a maximum greater than $0.2$
for $\frac{F}{\omega^{2}}=0.55$ (red line), then rapidly decreases to zero.
For $\frac{F}{\omega^{2}}=0.26$ (green line), the amplitude does not exceed
$0.05$ then it decreases to zero. It can be concluded that the potential and
the laser field parameters have a direct effect on all transmissions.

(a)

(b)

Figure 5: (Color online) Transmission probabilities as a function of potential
$V_{0}$ with $d=7$, $V_{1}=0$, $E=20$, $k_{y}=1$. (a): $\omega=1.4$, $F=0.9$
for total transmission $T$ (magenta line), the central transmission band
$T_{0}$ (blue line) and two first sidebands $T_{-1}$ (red line), $T_{1}$
(green line). (b): Transmission probability $T_{-1}$ with $F=1.8$,
$\omega=1.4$ (blue line), $\omega=1.8$ (red line) and $\omega=2.6$ (green
line).

The transmission probabilities for the first three bands are shown in Fig. 6
as a function of barrier width $d$, with $V_{0}=20$, $V_{1}=10$, $E=5$,
$k_{y}=1$, and $\omega=1.4$. In Fig. 6a, for $F=0.9$, we observe that all
transmissions oscillate in a sinusoidal manner with the same period but
different amplitudes, indicating a large difference between $T_{0}$ and
$T_{-1,1}$. The transmission $T_{0}$ is always more dominant, and it
oscillates in the vicinity of $1$. When we plotted the same transmissions for
$F=1.8$, we saw that their nature remained the same, but their amplitudes grew
larger without the spacing between $T_{0}$ and $T_{-1,1}$, as shown in Fig.
6b. The amplitude of $T_{0}$ decreases and that of $T_{-1,1}$ increases along
the $d$ axis, but the sum of the three transmissions, $T_{-1,0,1}$, would
never exceed unity, which indicates that the transmission in the central band
is important for short barriers. However, side band transmissions become
important for large barriers. We conclude that the laser parameters have an
effect on the amplitudes and the nature of the transmission probabilities.

(a)

(b)

Figure 6: (Color online) Transmission probabilities for the central band
$T_{0}$ (blue line) and two first sidebands $T_{-1}$ (red line), $T_{1}$
(green line) as a function of the barrier width $d$ with $V_{0}=20$,
$V_{1}=10$, $E=5$, $k_{y}=1$. $\omega=1.4$. (a): $F=0.8$ and (b): $F=1.8$.

##  V Conclusion

We theoretically investigated a system consisting of three graphene regions:
$1,2,3$, virgin graphene, and an intermediate region with a laser field
generated by an electric field of amplitude $F$ and frequency $\omega$ in the
presence of an inclined scalar potential generated between two identical
plates by two different potential generators $V_{0},V_{1}$. First, we solved
the eigenvalue equation to determine the wave functions corresponding to each
region. For this, we used the Floquet’s theory and the solution of Weber’s
differential equation. To determine the transmission probabilities, we
employed the boundary conditions and used transfer matrix approach, which give
an infinite order transfer matrix. To simplify, we limited our studies to the
first three bands, that are the central band $l=0$ and the two first side
bands $l=\pm 1$.

We plotted the transmission probabilities as a function of the physical
parameters $E$, $V_{0}$, $V_{1}$, $d$, and $\omega$. In summary, our numerical
results showed that the laser field plays an important role in the barrier.
The oscillations of this barrier over time generate several modes of
transmission probabilities. It is demonstrated that the reflection is complete
inside the barrier for all modes. It is found that the variation of the
frequency $\omega$ or the amplitude $F$ (i.e. $\frac{F}{\omega^{2}}$) changes
the form and the amplitude of transmissions. More importantly, we observed
that the transmission of the central band is always more dominant than the
sidebands, and the Klein zone is always present, especially for energies lower
than the barrier potential.

## References

  * [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * [2] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * [3] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * [4] S. Bhattacharjee and K. Sengupta, Phys. Rev. Lett. 97, 217001 (2006).
  * [5] J. S. Bunch, Y. Yaish, M. Brink, K. Bolotin, and P. L. McEuen, Nano Lett. 5, 2887 (2005).
  * [6] Claire Berger, Zhimin Song, Tianbo Li, Xuebin Li, Asmerom Y. Ogbazghi, Rui Feng, Zhenting Dai, Alexei N. Marchenkov, Edward H. Conrad, Phillip N. First, and Walt A. de Heer, J. Phys. Chem. B 108, 19912 (2004).
  * [7] Y. Zhang, Y. W. Tan, H. L. Stormer, and P. Kim, Nature 438, 201 (2005).
  * [8] S. V. Morozov, K. S. Novoselov, F. Schedin, D. Jiang, A. A. Firsov, and A. K. Geim, Phys. Rev. B 72, 201401 (2005).
  * [9] H. Haugen, D. Huertas-Hernando, A. Brataas, Phys. Rev. B 77, 115406 (2008).
  * [10] Z. H. Ni, T. Yu, Y. H. Lu, Y. Y. Wang, Y. P. Feng, and Z. X. Shen, ACS Nano 2, 2301 (2008).
  * [11] M. Y. Huang, H. G. Yan, C. Y. Chen, D. H. Song, T. F. Heinz, and J. Hone, Proc. Natl. Acad. Sci. U.S.A. 106, 7304 (2009).
  * [12] Y. Zahidi, I. Redouani, A. Jellal, and Hocine Bahlouli, Physica E 115, 113672 (2020).
  * [13] H. Bahlouli, E. B. Choubabi, A. El Mouhafid, and A. Jellal, Solid State Communications 151, 1309 (2011).
  * [14] M. Mekkaoui, A. Jellal, and H. Bahlouli, Solid State Communications 358, 114981 (2022).
  * [15] C. Sinha and R. Biswas, Appl. Phys. Lett. 100, 183107 (2012).
  * [16] R. Biswas and C. Sinha, Appl. Phys. 114, 183706 (2013).
  * [17] R. El Aitouni and A. Jellal, Phys. Lett. A. 447, 128288 (2022).
  * [18] H. L. Calvo, H. M. Pastawski, S. Roche, and L. E. F. Torres, Appl. Phys. Let. 98, 232103 (2011).
  * [19] S. E. Savel’ev and A. S. Alexandrov, Phys Rev. B. 84, 035428 (2011).
  * [20] J. T. Lu, F. H. Su, H. Wang, and X. H. Deng, Europhys. Lett. 95, 24003 (2011).
  * [21] M. Mekkaoui, R. El Kinani, and A. Jellal, Mater. Res. Express 6, 085013 (2019).
  * [22] C. W. J. Beenakker, Rev. Mod. Phys. 80, 1337 (2008).
  * [23] M. I. Katsnelson, K. S. Novoselov, A. K. Geim, Nat. Phys. 2, 620 (2006).
  * [24] R. Biswas, S. Maitty, and C. Sinha, Physica E 84, 235 (2016).
  * [25] E. Choubabi, A. Jellal, and M. Mekkaoui, Eur. Phys. J. B 92, 1 (2019).
  * [26] R. Loudon, The Quantum Theory of Light, 3rd ed. (Oxford University Press Inc, New York, 2000).
  * [27] Z. Gu, H. A. Fertig, D. P. Arovas, and A. Auerbach, Phys. Rev. Lett. 107, 216601 (2011).
  * [28] F. W. J. Olver, J. Res. Nat. Bur. Standards B 63, 131 (1959).
  * [29] A. Jellal, E. B. Choubabi, H. Bahlouli, and A. Aljaafari, J. Low Temp. Phys. 168, 40 (2012).


A Semblance-based Microseismic Event Detector for DAS Data–LABEL:lastpage

# A Semblance-based Microseismic Event Detector for DAS Data

Juan Porras1,2 Davide Pecci1,3 Gian Maria Bocchini4 Sonja Gaviano1 Michele De
Solda1 Katinka Tuinstra5 Federica Lanza5 Andrea Tognarelli1 Eusebio Stucchi1
and Francesco Grigoli1
1 Department of Earth Sciences  University of Pisa  _56124_ Pisa  Italy
2 Department of Earth Sciences  University of Geneva  _1205_ Geneva
Switzerland. Email<EMAIL_ADDRESS>
3 DESTeC  University of Pisa  _56124_ Pisa  Italy
4 Institute of Geology  Mineralogy and Geophysics  Ruhr University of Bochum
_44801_ Bochum  Germany
5 ETH-Zurich  Swiss Seismological Service (SED)  _CH8092_ Zurich  Switzerland

(Received 1998 December 18; in original form 1998 November 22; 1998)

###### keywords:

Distributed Acoustic Sensing, Microseismic Monitoring, Earthquake Detection.

††volume: 200

Distributed Acoustic Sensing (DAS) is becoming increasingly popular in
microseismic monitoring operations. This data acquisition technology converts
fiber-optic cables into dense arrays of seismic sensors that can sample the
seismic wavefield produced by active or passive sources with a high spatial
density, over distances ranging from a few hundred meters to tens of
kilometers. However, standard microseismic data analysis procedures have
several limitations when dealing with the high spatial (inter-sensor spacing
up to sub-meter scale) sampling rates of DAS systems. Here we propose a
semblance-based seismic event detection method that fully exploits the high
spatial sampling of the DAS data. The detector identifies seismic events by
computing waveform coherence of the seismic wavefield along geometrical
hyperbolic trajectories for different curvatures and positions of the vertex,
which are completely independent from external information (i.e. velocity
models). The method detects a seismic event when the coherence values overcome
a given threshold and satisfies our clustering criteria. We first validate our
method on synthetic data and then apply it to real data from the FORGE
geothermal experiment in Utah, USA. Our method detects about two times the
number of events obtained with a standard method when applied to 24h of data.

##  1 Introduction

Distributed Acoustic Sensing (DAS) is a fiber optic-based data acquisition
technology that is gaining popularity in a broad range of seismological
applications. From seismic monitoring in volcanic (Jousset et al., 2022),
volcano-glacial environments (Klaasen et al., 2021), offshore areas (Shinohara
et al., 2022) and aftershock scenarios (Li et al., 2021), (Zeng et al., 2022);
to downhole (Lellouch & Biondi, 2021), surface applications (Lindsey & Martin,
2021), traffic and railway monitoring (Wang et al., 2021). Different
industrial applications also benefit from the use of DAS, such as induced
seismicity monitoring (Webster et al. 2013; Karrenbach et al. 2019), CO2
sequestration monitoring (Daley et al., 2013), hydraulic stimulation
(Karrenbach et al., 2017) and geothermal sites (Lellouch et al. 2020; Lellouch
et al. 2021).

A DAS system operates by sending a laser pulse along a fiber-optic cable. Due
to small fluctuations (i.e. inhomogeneities) of the refractive index within
the fiber, a portion of the light pulse travelling through the cable is back-
scattered. The reflected pulse is sent back along the fiber to the source
device that also works as a receiver. This device is commonly known as
Interrogator Unit (IU). The IU exploits the Rayleigh back-scattering
principle, where the transmitted and the back-scattered pulse have the same
frequency, and the phase difference between the two can be translated to
dynamic strain along the fiber. When the fiber is unperturbed (i.e. no static
or dynamic deformation is occurring), the arrival phase of the back-scattered
light is constant. In the presence of an external perturbation (e.g. an
incident seismic wave), the fiber length changes and produces a corresponding
phase delay in the back-scattered pulse, which is recorded and converted by
the IU into strain or strain rate. More technical details about DAS
acquisition systems can be found in Lindsey et al. (2020) and Paitz et al.
(2021).

The use of DAS technology is constantly increasing in induced seismicity
monitoring with borehole installations (Lellouch & Biondi, 2021). In such
applications, the typical monitoring setting consists of a geophone chain
deployed in a deep borehole, allowing to detect a larger number of seismic
events with respect to conventional surface seismic networks (Maxwell et al.
2012; Rossi et al. 2023). Despite their extensive use, borehole arrays of
geophones have two main limitations: 1) their operational temperature may
strongly affect their installation depth, and 2) the need for a dedicated
monitoring well (or observation well) that increases operational costs. These
limitations are particularly pernicious in monitoring operations of induced
seismicity associated with Enhanced Geothermal Systems (EGS). In fact, the hot
dry rocks of an EGS may reach temperatures higher than 300 °C, hence deploying
conventional geophones down to the reservoir level is not feasible as the
electronics of these sensors is not designed to work in such high temperature
conditions (Zhidong et al., 2019). Last but not least, conventional
seismometer deployments sample the seismic wavefield in a sparse network of
point measurements, making these sensors of limited use to analyze the
complexity of the seismic wavefield in great detail. Alternatively, the fiber-
optic cable can resists to high temperatures (Zhidong et al., 2019). In
addition, DAS systems sample the seismic wavefield with a very high spatial
sampling (up to sub-meter scale) over tens of kilometers (Wang et al., 2020),
providing a more complete picture of the seismic wavefield. Thus, DAS
recordings can potentially enable us to detect features of the seismic
wavefields that would go unnoticed with conventional recordings. For example,
Lindsey et al. (2019) show that DAS recording were able to capture/image waves
from an on-shore earthquake reflected at previously unknown offshore faults
located below the deployed fiber-optic cable.

For these reasons, DAS is a promising technology for the monitoring of induced
seismicity associated to EGS operations. Another benefit of this technology is
that any production or stimulation well can be turned into an observation well
by just installing an optical-fiber behind the casing. This reduces the costs
of the monitoring infrastructure but also improves the earthquake detection
capability of DAS since the distance from the sensing fiber to the reservoir
is reduced. However, a downside of this technology is the massive amount of
produced data. An example takes place at the Frontier Observatory for Research
in Geothermal Energy (FORGE) site in Utah, USA. This is a dedicated
underground field laboratory aiming at developing, testing, and accelerating
breakthroughs in EGS technologies to advance the uptake of geothermal
resources around the world (McLennan, 2022) (Fig. 1). At FORGE there have been
several DAS acquisition campaigns. For instance, the one in April 2022
(McLennan, 2022), which provides the real dataset for this study, used two
observation wells, namely 78A-32, and 78B-32. The interrogated sensing
optical-fiber was 2482 m long. DAS data were acquired using the following
settings: a gauge length of 10 m, a channel spacing of 1 m, and laser launch
rate of 20 kHz. Data were saved with a sampling frequency of 4 kHz and in 15
second chunks. Such a configuration produced about 1.3 TB of data per day.

Figure 1: Location map of the FORGE site in Utah, USA. The triangles represent
the different wells used for the third stage of the April 2022 stimulation
campaign. Well 16A(78)-32 corresponds to the injection well, 58-32 (located at
905 m from the injection well) is a monitoring well with a geophone array used
by Geo Energy Suisse to build a seismic catalog. Wells 78A-32 and 78B-32
(located at 1308 m from the injection well) host the fiber-optic of the DAS
system. The data used in this study was recorded at well 78B-32. The well-
scheme illustrates the design of the DAS system used to monitor the
stimulation campaign during April 2022. Wells 78A-32 and 78B-32 are separated
by a surface distance of 85 m and connected to allow for simultaneous data
acquisition using a single DAS interrogator unit (IU). The system starts
recording at the bottom of well 78A-32 and the end of the array is at 1193 m
below the surface in well 78B-32.

The large datasets produced by DAS systems highlights the problem of data
storage. This is particularly true when dealing with weeks or months long
seismic monitoring campaigns. A potential solution that could partially solve
this challenging problem is to store only the target data, such as waveforms
containing only the signal of the seismic events of interest. Unfortunately,
standard pick-based seismological techniques do not exploit the high spatial
density of DAS data and current seismic event detectors may miss a non
negligible number of seismic events. For this reason, using standard detectors
during seismic monitoring campaigns with DAS, increases the risk of
permanently losing useful information. Another downside of DAS applications is
its single-component nature. That is, a DAS system measures only along the
axial direction of the fiber and, in borehole geometries, it provides no
azimuthal information for earthquake recordings (Tuinstra et al., 2023).

Since the use of DAS technology for earthquake monitoring is a yet-to-be-
established practice, few algorithms for earthquake detection have been
developed for DAS data. For instance, Lellouch et al. (2019) developed a
detection method based on a single-parameter scan of the incidence angle and
measurement of data coherence along different possible travel-time curves,
given a seismic velocity model. Seismic detection methods based on the
measurement of local data coherence have been proposed also for large-N dense
networks (Li et al., 2018). Other methodologies of earthquake detection for
DAS data based on Convolutional Neural Networks (CNN) have been developed
(Binder & Chakraborty 2019; Stork et al.2020, Huot et al. 2022 and Zhu et al.
2023). In our opinion, a seismic event detector suitable for DAS data must
meet the following criteria: 1) It must exploit the high spatial density of
DAS, 2) it must be sensitive but robust at the same time (i.e. low missed and
low false event detection rates) and 3) it needs to be computationally fast to
be applied in real-time or near real-time.

In this work we provide a solution to this important problem by implementing
HECTOR (coHerence-based Earthquake deteCTOR), a semblance-based seismic event
detection method that fully exploits the characteristics of DAS data. We first
introduce the theoretical background at the base of our earthquake detector.
Next, to evaluate the performance of HECTOR, we test it both on synthetic and
real DAS microseismic data acquired during the 2022 EGS stimulation campaign
at the Utah FORGE site for which a DAS-based seismic catalogue has been
already built by Silixa (LLC, 2022), this is our reference catalogue in this
study. Finally, we discuss the performance of the algorithm by comparing our
results with the reference catalogue.

##  2 Data and Methods

###  2.1 Detection method

The detection methodology consists of two parts: 1) the first one evaluates
the spatial coherence of the seismic waveforms using the semblance function;
and 2) the second part performs a clustering analysis of the coherence time
series obtained from the first part of the methodology to identify
microseismic events.

####  2.1.1 Waveform Coherence

HECTOR exploits the high spatial sampling of DAS for real-time microseismic
monitoring. In particular, for a linear segment of the optical-fiber, HECTOR
evaluates the coherence of the seismic waveforms along hyperbolic trajectories
by using the Semblance function (Neidell & Taner, 1971). The equation of the
hyperbola we use to calculate the coherence is defined as:

$t_{i}(T,X,C)^{2}=T^{2}+\frac{(x_{i}-X)^{2}}{C^{2}}$ (1)

where, $t$ is the time of the trace recorded at the position $x_{i}$ (along
the fiber axis), while the coefficients $X$, $T$ and $C$ represent
respectively the spatial offset with respect to the vertex of the hyperbola,
the time offset, and the curvature. Since we are interested only on positive
values, the previous equation can be written as:

$t_{i}(T,X,C)=\sqrt{T^{2}+\frac{(x_{i}-X)^{2}}{C^{2}}}$ (2)

Finally, the Semblance function $S$ can be written as:

$S(T,X,C)=\frac{\sum_{j=1}^{N}\left(\sum_{i=1}^{M}A(t_{ij})\right)^{2}}{M\sum_{j=1}^{N}\sum_{i=1}^{M}A(t_{ij})^{2}}\
\ \ \mathrm{with}\ \ \ t_{ij}=t_{i}(T,X,C)+j\mathrm{d}t$ (3)

Where $M$ is the total number of seismic traces (i.e. the number of the
recording elements of the fiber), $N$ is the length of the sample window, and
$A$ is the amplitude for the $i^{th}$ trace at the time $t_{ij}$. Depending on
the application, we can choose to use either raw or processed waveforms (e.g.
envelope) to calculate the semblance. The Semblance value ranges between 0 and
1, where 1 indicates perfect coherence and 0 no coherence of the signals from
the different traces. The detector identifies seismic events by evaluating the
waveform coherence along geometrical hyperbolic windows (sliding window
hereafter, Fig. 2.c) while changing their curvature and the position of the
vertex over the data space.

Figure 2: Sketch of the semblance-based detection method. a) represents the
coherence time series, which overcomes the detection threshold in the presence
of energy signal. b) shows a volume of 2D semblance matrices resulting from
the scanning of the waveform coherence along geometrical hyperbolic shapes.
Each semblance panel results from the scanning at each specific X position
along the spatial offset. c) illustrates the range of geometrical hyperbolas
(coefficients $C_{1}$ to $C_{N}$) at each $X$ and $T$ used to measure the
waveform coherence along a finite-width data window ($t_{j}$ to $t_{j}+\Delta
t$) (from now on sliding window) illustrated by the orange hyperbolic area.
The coefficients $X$ and $T$ allow for the scanning of events along the
spatial and temporal axes, respectively. In the absence of an energy signal,
the method returns low coherence values. However, in the presence of a
hyperbolic event, it is detected as a high coherence region, which, if
overcomes the detection threshold and satisfy our clustering criteria, is
declared as a seismic event.

The algorithm generates a volume of 2D semblance matrices (Fig. 2.b) whose
number of matrices and the rows and columns of each matrix is given by the
number of $X$, $C$ and $T$ coefficients, respectively. The latter is decided
upon the required resolution and processing times. Finally, the selected
semblance matrix is the one with the highest coherence values.

The hyperbolic trajectories that do not follow the seismic wavefield will have
an overall semblance value close to zero since it is equivalent to the sum of
random noise. On the contrary, the hyperbolic trajectories that do follow the
seismic wavefield, regardless the wave phase, will tend to have higher
coherence values. The hyperbolic trajectories used to measure the waveform
coherence are completely independent from external information (i.e. velocity
model and hypocentral location), therefore, our method uses the waveform data
as the only input. The output is a time series of coherence values that is the
result of squaring and summing the coherence values along the columns of the
selected semblance matrix. Squaring the values of the semblance matrix allows
us to reduce even more the small coherence values related to the random noise
while increasing the difference with higher coherence values, potentially
associated with seismic events. Each column represents a time step in the
computation of the semblance matrix and the number of columns depends on the
number of $T$ coefficients (Fig. 2.a). The coherence time series is a 1D
expression of the semblance matrix, obtained in the first part of the
detection methodology. This outcome is used as input for the second part of
the detection methodology, where we perform a clustering of the coherence
values that exceed a detection threshold to identify potential detections of
microseismic events.

It is important to note that applying the Semblance function directly to raw
seismic waveforms (i.e. zero mean traces) can lead, in some cases, to spurious
results. This is due to the non-isotropic radiation pattern of an earthquake
source. In such cases, waveforms recorded along the fiber may have reversed
polarities, thus reducing the overall coherence produced by the stacking
process. To avoid this issue one may use normalized STA/LTA functions based on
energy (a non-negative function), that mitigates the effect of the radiation
pattern (Grigoli et al., 2013). However, stacking non-negative functions (such
as the energy or the envelope of each trace) reduces the radiation pattern’s
impact but does not suppress the random noise. Random noise suppression can be
achieved by zero-mean stacking functions , such as raw waveforms. In the FORGE
case, the length of the fiber is limited compared to the distance from the
stimulation site. Such a condition allows the successful use of raw waveforms
without further data processing and without having the problem of the
radiation pattern effect (due to the limited sampling of the focal sphere).

####  2.1.2 Clustering Analysis

We calculate the detection threshold as the mean of the coherence time series
after removing the 5% of the lowest and largest values. The clustering
approach requires the calibration of two parameters: 1) the minimum number of
samples required to form a cluster, and 2) the maximum number of adjacent
samples below the detection threshold to consider a single cluster. Clusters
that meet the previous parameters are considered as potential microseismic
events where the detection time corresponds to the time of the first sample in
the cluster (Fig. 3).

To refine the initial list of detections from the cluster analysis, we apply a
signal-to-noise-ratio (SNR) criterion and calculate the RMS of the signal and
noise windows as follows:

$E_{RMS}=\sqrt{\frac{\sum_{i=1}^{N}w_{i}^{2}}{N}}$ (4)

Where $N$ is the length of the time window in samples and $w_{i}$ is the
$i$-th sample of the waveform. Then, the SNR is computed as:

$SNR=20log10\frac{E_{RMS{signal}}}{E_{RMS{noise}}}$ (5)

The SNR threshold used to refine the initial list of detections will depend on
the data quality and noise levels and can be manually tuned to enhance the
number of real detections and minimize the number of false detections.

Figure 3: Schematic representation of the clustering approach used to
individuate seismic events from the coherence time series (Fig. 2). A seismic
event is returned when a cluster contains a minimum number of samples above
the detection threshold and passes a SNR control (see text). The number of
samples within the coherence time series is equal to the number of columns of
the semblance matrix. The value of each sample of the coherence time series is
equal to the sum along the columns of the coherence matrix of the squared
coherence values (Fig. 2). When an event is detected, the detection time of
the event corresponds to the time of the first sample in the coherence time
series. When the cluster size is larger than a time longer than the maximum
duration of the targeted microseismicity (MDTM), we inspect the possibility
that a given cluster contains two events (see sub-section Clustering
Analysis).

The signal and the noise windows are represented by a certain number of
coherence samples after and before a detection time, respectively. The noise
window is shorter than the signal window. We retain the detection as an event
if the SNR is larger than a certain threshold, otherwise we reject it. In the
case of consecutive coherence samples exceeding the detection threshold for a
time longer than the maximum duration expected for the targeted
microseismicity, we inspect the possibility of the presence of two
microseismic events within the same cluster of coherence samples. Multiple
microseismic events falling within the same cluster are expected when the
inter-event distances are very short (Fig. 3). In the case of clusters of
coherence samples exceeding the duration of targeted microseismic events, we
search for the largest positive step (i.e. the largest difference in
subsequent coherence values) in the coherence values and consider it as a
potential microseismic event. To confirm the event, we apply a SNR criterion
like the one applied for the first detection. The method is currently unable
to detect more than two events should they fall within the same cluster of the
coherence time series.

It is worth mentioning that this methodology is developed specifically to
detect microseismic events monitored in straight fibers where the seismic
wavefield can be approximated as hyperbolic events. However, HECTOR also
allows to detect events originating in the far field if low curvature
parameters are provided to match plane waves and a proper tuning of the
clustering parameters.

In the next sections, we test and validate the performance of our detector
(HECTOR) by applying it both to synthetic and real data.

###  2.2 Data pre-processing

A sequence of pre-processing steps are applied to the data before running the
detector. In this work we defined an automatic denoising workflow (Fig. 4)
that is crucial to improve the capability of our detector especially for weak
or very low-magnitude events.

The denoising workflow consists of the following steps: 1) Application of a
low-pass filter (before downsampling to avoid aliasing effects) and
downsampling of the traces, the latter to enhance computational speed. 2)
Removal of the linear and mean trends of each trace of the array. 3)
Normalization of the trace amplitudes by its corresponding maximum to reduce
the effect of the geometrical spreading, receiver coupling and nonlinear
effects (see Miah & Potter., 2017 for a detailed review). 4) Application of a
band-pass filter to remove the high frequency and low frequency noise on the
traces and isolate the frequency band of interest for the events in the
catalogue. 5) Application of a Frequency-Wavenumber (FK) filter to attenuate
any coherent noise. The FORGE dataset is characterized by strong coherent
noise parallel to the fiber. The FK filter is an optimal array seismology
technique used to attenuate this kind of unwanted energy.

Fig. 4 shows an example of a low-amplitude microseismic event (maximum
amplitude 44 $n\epsilon/s$) before and after applying the denoising workflow.
The event is slightly visible in the time domain and in the FK spectrum of the
raw data (Fig. 4 a and b). However, after applying the denoising workflow, the
coherent noise along the $k=0$ wavenumber is strongly attenuated as well as
the high frequency random noise, highlighting the spatial coherence of the
seismic event (Fig. 4 c and d).

Figure 4: Example of the application of the denoising workflow. a) Raw DAS
data with strong coherent noise parallel to the fiber-optic cable. b) FK
spectrum of the raw DAS data. The strong coherent noise is mapped along the
wavenumber zero and hinders the spectrum of the energy signal. c) Denoised DAS
data. d) FK spectrum of the denoised DAS data. The color scale of the FK
spectra are normalized to 1 but saturated to 0.5 for visualization purposes.

###  2.3 Data

We first generated a synthetic DAS dataset of 36 events sampled at 500 Hz,
with an inter-time distance of 5 seconds, central frequency of 50 Hz and local
Magnitudes (Ml) ranging between -2.0 to 1.5 by using the software Salvus
(Afanasiev et al., 2019) and the velocity model available in Lellouch et al.
(2020). The synthetic dataset resembles real DAS data recorded at FORGE, and
allows us to have full control over the earthquake source parameters, test the
performance of HECTOR and validate our methodology.

To build the dataset we first generated the synthetic waveforms (Fig. 5) for a
reference event recorded during the stimulation campaign in April 2019. Using
the earthquake catalogue available in Lellouch et al. (2020), we defined as
reference an event of Ml 0.91 that occurred on April 23th, 2019 at 21:32:09
UTC, at hypocentral distance of 3390 m from the bottom of the fiber. We
simulated the synthetic strain waveforms for a fiber cable of 1.2 km with 1 m
inter-channel spacing (as in the real dataset). We then converted the
waveforms in strain-rate and corrected the amplitudes by scaling them with the
waveform amplitudes of the corresponding real event. The scaling of the
amplitudes is necessary to achieve synthetics that are closer to the recorded
signals.

Finally, to simulate the full set of events with different magnitudes we re-
scaled the amplitudes of the synthetic waveforms using the following relation:

$A_{i}^{M}(t)=A_{i}^{R}(t)\frac{10^{M}}{10^{M^{R}}}$ (6)

where $A_{i}^{R}(t)$ are the waveform amplitudes of the reference event with
magnitude $M^{R}$ (here $M^{R}$ = 0.91), while $A_{i}^{M}(t)$ are the scaled
amplitudes for an event of magnitude $M$. The index $i$ refers to the $i$th
DAS sensor. Finally, to obtain a dataset as similar as possible to the real
data, we superimposed Gaussian distributed noise to the synthetic traces. We
first calculate the real noise as the mean of the RMS of four 15 seconds
chunks of data (i.e. four files). We then calculate the standard deviation of
the data in the four 15 seconds chunks to simulate Gaussian distributed noise.
Finally, we apply to the synthetic traces the denoising workflow of section
2.2 to measure the maximum amplitudes before calculating the noise. This
synthetic dataset aims to reproduce different SNR conditions rather than
reproducing the exact amplitudes of events of different magnitudes which is
beyond the scopes of this study.

We then evaluated the performance of HECTOR by applying it to real data
acquired at FORGE during the stimulation campaign in April 2022. We used data
acquired with the Silixa’s Carina sensing system along two simultaneously
interrogated observation wells, namely 78A-32, and 78B-32. As mentioned
before, the interrogated optical-fiber had a length of 2482 m and DAS data
were acquired using a gauge length of 10 m and a channel spacing of 1 m. Data
were saved with a sampling frequency of 4 kHz and in files of 15 second
chunks. Additional details on the DAS acquisition can be found in the
Stimulation Silixa Microseismic Report (LLC, 2022).

Figure 5: a) 15 seconds of real DAS data containing the Ml 0.91 earthquake
occurred on April 23th, 2019 at 21:32:09 UTC, located at a hypocentral
distance of 3390 m from the very end of the fiber, used to build the synthetic
waveforms. b) 15 seconds of synthetic DAS data. The synthetic event has an
amplitude equivalent to a Ml 0.91 event and it is located at a depth of 1600 m
and epicentral distance of 3000 m from the well.

We focus on data acquired during the third stimulation stage, starting on
April 21st 2022 at 13:48:22 UTC, that generated the largest number of
microseismic events among the three stimulation stages performed in April
2022. We downloaded 24 hours of DAS data (1.3 TB), from April 21st 2022 at
13:00 UTC to April 22st 2022 at 13:00 UTC. Silixa’s catalogue contains 1199
events within the investigated time interval (LLC, 2022).

##  3 Results

In this section we describe the parameters chosen for each test and report the
detection results for the synthetic and real datasets, respectively.

###  3.1 Results for the synthetic dataset

A required step in the development of novel seismological algorithms consists
in extensive testing with synthetic data.

Here, for our synthetic tests, we apply HECTOR with the following parameters:
sliding windows ($N$) of 20 samples and steps ($T$) of 10, and hyperbolic
curvature values ($C_{1}$ to $C_{N}$) between 50 and 200 with steps of 5. The
curvature parameters are dependent on the distance of the events from the
fiber. The size and steps of the sliding window are decided upon the sampling
rate of the data and the required resolution. We fixed the lateral search
parameter ($X$) to be the depth of the deepest sensor which is a good
assumption when the microseismic events are deeper than the bottom of the
fiber. We set to eight the minimum number of samples over the detection
threshold to declare an event and to three the maximum number of consecutive
samples below the detection threshold to identify the synthetic microseismic
events. The previous parameters were defined based on a trial and error
procedure.

Figure 6: 60 seconds of synthetic DAS data containing events from Ml -0.8 to
0.3 (from right to left) with overlapped Gaussian noise. a) Coherence time
series where the starting time of each event is indicated by the black
vertical dashed-line. b) Semblance matrix resulting from the scanning of the
waveform coherence along geometrical hyperbolic shapes. c) Synthetic event
waveforms separated by 5 second windows. The dashed black line indicates the
channel shown in panel d. d) Waveform from a single channel that contains the
events from Ml -0.8 to 0.3 (from right to left). Top-right of the event
waveforms is reported the SNR calculated on the shown trace. SNR is estimated
using a noise window of 1 second before the P-wave arrival and a signal window
of 1.5 seconds after the P-wave arrival.

We report the results of the synthetic test in Fig. 6 and Fig. S1 of the
Supplementary material. The application to the synthetic dataset clearly shows
that HECTOR can successfully detect microseismic events with very low SNR or
even microseismic events hidden in the noise if the signal is coherent across
multiple channels of the fiber (Fig. 6). For the synthetic data we also
observe a good agreement between the P-wave arrival of the events (Fig. 6c)
and the detection time (Fig. 6a). Microseismic events with lower amplitudes
(Fig. 6cd) result in lower coherence values (Fig. 6ab). However, we are unable
to detect signals with extremely low amplitudes with respect to the noise
amplitude as shown in Fig. S1 of the supplement. After successfully validating
HECTOR with the synthetic dataset we apply it to real data.

###  3.2 Results for the real dataset

Figure 7:  a) Coherence time series for 3 events in a 15 seconds window.
Detection times are: 17:06:28.5, 17:06:29.9 and 17:06:34.8 UTC on April 21st,
2022. Dashed vertical black lines indicate the detection time of the
microseismic events. b) Semblance matrix in which the three regions of maxima
correspond to the three microseismic events detected by HECTOR c) Real DAS
data from FORGE (well 78B-32) that contains three microseismic events with
varying amplitudes. The first two events occur between 4 and 6 seconds while
the third event, with the smallest amplitude among the three, occur between 10
and 11 seconds. DAS data in panel $c$ is already filtered following the
denoising workflow (section 2.2). Figure 8: Example of two microseismic events
present in the catalog from this study (detected by HECTOR) and not in
Silixa’s catalog. a) Event with detection time on April 21st 2022, at
14:46:11.3 UTC. b) Event with detection time on April 21st 2022, at 17:47:38.3
UTC. These events have not been identified by the STA/LTA algorithm used by
Silixa. Figure 9: Detections from HECTOR applied to 24h of real data from
FORGE and comparisons with an existing catalogue. Histogram bins are of 30
minutes. a) Common detections between this study and the Silixa’s catalogue
(gray), detections obtained only from this study (red), and the detections
present only in the Silixa’s catalogue (blue). b) Total detections from this
study (blue), total detections in the Silixa’s catalogue (orange), total false
detections from this study (gray), and total false detections from Silixa’s
catalogue (green).

In order to apply the detector to the real data, we calibrated the detection
and clustering parameters on a set of 50 files containing: single events,
multiple events with small inter-distances (Fig. 7), or only noise (Fig. S2).
We observe that the typical duration of microseismic events recorded along the
fiber-optic cable is 0.6-0.8 seconds. To enhance computational speed, we used
only the channel interval from 1350 to 2384 (Fig. 7), corresponding to the
optical-fiber cable located along the observational well 78B-32. We omit the
topmost channels from the well 78B-32 (from 1215 to 1350) due to the high
noise levels near the surface. We downsampled the real data to 600 Hz, applied
a band-pass filter between 10 to 250 Hz and a FK filter to attenuate coherent
linear noise and unwanted energy.

The calibrated detection parameters for the evaluation of the waveform
coherence consist of a sliding window of 20 samples with a time search step
($T$) of 10 samples and the curvature values ($C_{1}$ to $C_{N}$) from 100 to
180 with a step of 5. The lateral search was fixed to match the depth of the
deepest sensor as in the case of the synthetics. This assumption is valid in
the specific case of FORGE where the events are deeper than the bottom of the
fiber.

Regarding the clustering parameters for the detection of microseismic events,
we set to 10 the minimum number of samples above the detection threshold to
identify a cluster (i.e. potential detection). If two or more consecutive
clusters (or samples above the threshold) are separated by less than two
samples below the detection threshold, they are grouped into the same cluster.
As mentioned in the methodology, we apply a SNR criterion to refine the list
of detections. We defined a signal window to be 30 samples (corresponding to
0.5 seconds) of the coherence time series and the noise window to be 20
samples (from 2 to 22 coherence samples before the detection time,
corresponding to about 0.3 s) (Fig. 3). The second pair of signal and noise
windows to explore the presence of two events within the same cluster are
equal to the first 20 samples of the coherence time series after the
additional detection time, while the noise window includes 15 samples, namely
from 2 to 17 samples before the additional detection time, Fig. 3). We kept
the detection as an event if the SNR was larger than 4.

After the calibration, we processed the 24 hours of real data and initially
detected 2236 events (before de-clustering). As comparison, Silixa’s catalogue
contains 1199 detections for the same time interval (before de-clustering).
Their catalogue contains microseismic events detected using an STA/LTA method
applied to the fiber-optic cables deployed in both the observation wells.

To better show the performance of our algorithm, we made a rigorous comparison
of the events in our catalogue with those in the Silixa’s catalogue (LLC,
2022). We removed possible duplicated events by de-clustering both catalogues
with the removal of the second detection when two detections are less or equal
than 0.7 seconds apart. We retained 2187 and 1169 events, respectively. Then,
we searched for common events in the catalogues by considering a time interval
of ±0.6 s from the detection time of the events in our catalogue. We found
1086 common events, while 83 events are only contained in the Silixa’s
catalogue, and 1101 events are only contained in our catalogue (Fig. 8). We
manually inspected all the detected events and classified as ”false detection”
6 out of the 1086 events common to both catalogues, 64 out of the 83 events
only present in the Silixa’s catalogue, and 116 of the 1101 events only
present in our catalogue. We classified as ”false detection” the coda of
microseismic events, events with not clearly visible P and/or S arrivals,
false detections, and distant seismic events (Fig. S3).

In summary, our catalogue contains 2187 microseismic events (Fig. 9), from
which 122 are classified as ”false detection” (5.6 %), while the Silixa’s
catalogue contains 1169 microseismic events with 70 of them classified as
”false detection” (6.0 %). The reported total number of events also includes
common detections to both catalogues. The pre-processing workflow and the
semblance-based detection method proposed in this study enable the detection
of a number of events that is significantly larger (almost double) than those
detected with a standard STA/LTA method.

##  4 Discussion and Conclusions

In this work we developed HECTOR, a microseismic event detection method based
on the analysis of the spatial coherence of seismic wavefields along
geometrical hyperbolic trajectories.

After defining a proper DAS denoising workflow to improve the detection rate,
we first validated our detection method with synthetic DAS waveforms that
resemble the real data of the FORGE geothermal site recorded during the April
2019 campaign, and established a detection rate of microseismic events close
to a SNR of 0.1 (Figs. 5 and S4). We then applied the detector to 24 hours of
real data acquired during the third stimulation stage at FORGE during the
April 2022 acquisition campaign. With our method, we detected a total of 2187
microseismic events, from which only 5.6 % are false detections (2065 real
detections).

We have demonstrated that our detection method outpaces traditional pick-based
detection methods like STA/LTA for microseismic monitoring with DAS data, as
we almost doubled the number of detected events using the same dataset. During
the same period, Geo Energie Suisse, using borehole geophones deployed in well
58-32, closer to the stimulation well (905 m) with respect to well 78B-32
(1308 m) used in this study (Fig. 1), derived a catalogue of 1431 well-
locatable events with magnitudes ranging from 0.52 to -1.8 (Dyer et al.,
2022a). It is worth mentioning that Geo Energie Suisse reported a preliminary
total number of about 18400 detections during the stimulation stage 3 (Dyer et
al., 2022b). A direct comparison with the seismic catalogues obtained with
borehole geophones, as we did for the catalogue obtained from DAS data ((LLC,
2022)), is challenging and beyond the scope of this paper because of the
different locations of the monitoring wells.

In this study we provide a methodology that exploits the complexity of the
seismic wavefield in such a detail only possible to be recorded with DAS, to
perform an accurate detection of microseismicity, often characterized by a
massive number of small, noise contaminated, seismic events with short inter-
event times and without the need of external information (i.e. velocity
models). Our approach supports the potential benefits of turning any
production or stimulation well into an observation well too by just installing
an optical-fiber behind the casing, which would significantly reduce
operational costs.

Our detection algorithm has been developed specifically to detect microseismic
events monitored in straight fibers, where the seismic wavefield can be
approximated as hyperbolic events. However, HECTOR also allows to detect
events originating in the far field if low curvature parameters are provided
to match plane waves and a proper tuning of the clustering parameters.

The parameter $X$ allows for the scanning of microseismic events along
hyperbolic trajectories at different positions along the spatial offset
(parallel to the fiber-optic cable), even beyond the distance covered by the
DAS system. In the cases where an $X$ position, which represents the vertex of
the hyperbola, is beyond the fiber-optic cable, the hyperbolic trajectory used
to measure the spatial coherence will be only a part of the hyperbola (i.e.
the part covered by the fiber-optic cable). For each $X$ value, a semblance
matrix ($C-T$) is computed. Hence, we compute $N$ semblance matrices, where
$N$ represents the number of $X$ values considered in the lateral search
parameter. Among these $C-T$ semblance matrices, it is then selected the one
having the highest maximum coherence value. In this way the detector can be
applied to approximately straight fibers independently on the hypocentral
position of the event. However, lateral search also implies calculating a
number of $C-T$ semblance matrices equal to all the possible vertex positions
along the spatial offset, resulting in a higher computational cost. This
limitation will be improved by introducing parallelization in a future version
of the algorithm. Nonetheless, with the parameters defined in this study, the
algorithm requires a computational time of approximately 0.25 seconds per each
second of data, measured in a laptop with an Intel quad-core i7 processor and
16 GB of Random Access Memory (RAM).

Fig. S4 illustrates the effect of the parameter $X$ for three synthetic events
at different hypocentral locations with respect to a vertical well of 1 km
depth. We see that in all the cases, the events are successfully detected. The
event in Fig. S4c highlights how the detector successfully identifies both P
and S arrivals. Besides them, it also shows at least another detection related
to a best-fitting hyperbola having a lower curvature coefficient. This further
detection can be explained by considering that the algorithm tests several
potential hyperbolas having different curvatures. For a homogeneous and
isotropic half-space having constant velocity, the signal at the fiber would
be a perfect hyperbola. However, when dealing with 1D velocity models, the
signal recorded at the fiber can be seen as a combination of different
hyperbolas with different curvatures. As a result, for the same event and
seismic phase, there may be more than one detection related to hyperbolas
having different curvatures and fitting different portions of the signal.
Hence, multiple detections associated with the same event can be filtered out
by imposing a minimum interval between consecutive detections. As we observe
in Figure S5a,b, sometimes the detection is a few hundred milliseconds
anticipated. There may be two main reasons for this time difference. The first
is due to the best-fitting hyperbola’s vertex position along the time axis.
Given a borehole installation and a source beneath the fiber (as in Fig. S5b),
the best-fitting hyperbola vertex is located at a depth D greater than the
fiber depth and at a time instance that precedes the arrival time of the
signal at the fiber. The detection time corresponds to the hyperbola vertex
position along the time-axis. Hence, the time difference between the vertex
time position and the signal arrival time at the fiber is the observed delay.
On the other hand, in Figure S5a, the anticipated detection is due to the
clustering algorithm. The detector identifies clusters of consecutive
coherence samples whose amplitude exceeds a threshold. If a cluster is
classified as an event, the detection is set at the cluster starting time,
which may not correspond to the maximum of the coherence times series for that
cluster.

Despite the aforementioned limitations, we have demonstrated that our
waveform-based and clustering detection method outpaces traditional pick-based
detection methods and makes of HECTOR a suitable methodology for DAS-based
real-time microseismicity monitoring in industrial operations.

###### Acknowledgements.

We would like to thank the FORGE project for making the DAS data used in this
study freely available (McLennan, 2022). This work is part of a collaboration
between the Department of Earth Sciences of the University of Pisa and the
Swiss Seismological Service in the framework of the De-Risking Enhanced
Geothermal Energy Projects.(Innovation for DEEPs). DEEP is subsidized through
the Cofund GEOTHERMICA, which is supported by the European Union’s HORIZON
2020 programme for research, technological development, and demonstration
under Grant Agreement Number 731117. This work is also supported by the
University of Pisa under the Progetto di Ricerca di Ateneo PRA Fluid migration
in the upper crust: from natural hazards to geo-resources (project number 274
PRA 2022 66). G.M.B was funded by the Volkswagen Foundation and by the RDSME
of the Ruhr University of Bochum.

The list of detections from this study and the comparison with the Silixa’s
detections for the same period is freely available in (Bocchini et al., 2023).
The detector will be made publicly available on GitHub upon acceptance of the
paper.

## References

  * Afanasiev et al. (2019) Afanasiev, M., Boehm, C., van Driel, M., Krischer, L., Rietmann, M., May, D. A., Knepley, M. G., & Fichtner, A., 2019.  Modular and flexible spectral-element waveform modelling in two and three dimensions, Geophysical Journal International, 216(3), 1675–1692.
  * Binder & Chakraborty (2019) Binder, G. & Chakraborty, D., 2019.  Detecting microseismic events in downhole distributed acoustic sensing data using convolutional neural networks, in SEG Technical Program Expanded Abstracts 2019, pp. 4864–4868, Society of Exploration Geophysicists.
  * Bocchini et al. (2023) Bocchini, G. M., Porras, J., Pecci, D., & Grigoli, F., 2023.  Microseismic event catalogue for the april 2022 third stimulation at forge enhanced geothermal system (utah, usa) generated applying a semblance-based detection method to das data.
  * Daley et al. (2013) Daley, T. M., Freifeld, B. M., Ajo-Franklin, J., Dou, S., Pevzner, R., Shulakova, V., Kashikar, S., Miller, D. E., Goetz, J., Henninges, J., et al., 2013\.  Field testing of fiber-optic distributed acoustic sensing (das) for subsurface seismic monitoring, The Leading Edge, 32(6), 699–706.
  * Dyer et al. (2022a) Dyer, B., Karvounis, D., & Bethmann, F., 2022a.  Seismic data from the well 16a(78)-32 stimulation april, 2022.
  * Dyer et al. (2022b) Dyer, B., Karvounis, D., & Bethmann, F., 2022b.  Utah forge seismic events related to the april, 2022 well 16a(78)-32 stimulation.
  * Grigoli et al. (2013) Grigoli, F., Cesca, S., Vassallo, M., & Dahm, T., 2013.  Automated seismic event location by travel-time stacking: An application to mining induced seismicity, Seismological Research Letters, 84(4), 666–677.
  * Huot et al. (2022) Huot, F., Lellouch, A., Given, P., Luo, B., Clapp, R. G., Nemeth, T., Nihei, K. T., & Biondi, B. L., 2022.  Detection and characterization of microseismic events from fiber-optic das data using deep learning, Seismological Society of America, 93(5), 2543–2553.
  * Jousset et al. (2022) Jousset, P., Currenti, G., Schwarz, B., Chalari, A., Tilmann, F., Reinsch, T., Zuccarello, L., Privitera, E., & Krawczyk, C. M., 2022.  Fibre optic distributed acoustic sensing of volcanic events, Nature communications, 13(1), 1753.
  * Karrenbach et al. (2017) Karrenbach, M., Kahn, D., Cole, S., Ridge, A., Boone, K., Rich, J., Silver, K., & Langton, D., 2017.  Hydraulic-fracturing-induced strain and microseismic using in situ distributed fiber-optic sensing, The Leading Edge, 36(10), 837–844.
  * Karrenbach et al. (2019) Karrenbach, M., Cole, S., Ridge, A., Boone, K., Kahn, D., Rich, J., Silver, K., & Langton, D., 2019.  Fiber-optic distributed acoustic sensing of microseismicity, strain and temperature during hydraulic fracturing, Geophysics, 84(1), D11–D23.
  * Klaasen et al. (2021) Klaasen, S., Paitz, P., Lindner, N., Dettmer, J., & Fichtner, A., 2021.  Distributed acoustic sensing in volcano-glacial environments—mount meager, british columbia, Journal of Geophysical Research: Solid Earth, 126(11), e2021JB022358.
  * Lellouch & Biondi (2021) Lellouch, A. & Biondi, B. L., 2021.  Seismic applications of downhole das, Sensors, 21(9), 2897\.
  * Lellouch et al. (2019) Lellouch, A., Yuan, S., Ellsworth, W. L., & Biondi, B., 2019.  Velocity-based earthquake detection using downhole distributed acoustic sensing—examples from the san andreas fault observatory at depthvelocity-based earthquake detection using downhole distributed acoustic sensing, Bulletin of the Seismological Society of America, 109(6), 2491–2500.
  * Lellouch et al. (2020) Lellouch, A., Lindsey, N. J., Ellsworth, W. L., & Biondi, B. L., 2020.  Comparison between distributed acoustic sensing and geophones: Downhole microseismic monitoring of the forge geothermal experiment, Seismological Society of America, 91(6), 3256–3268.
  * Lellouch et al. (2021) Lellouch, A., Schultz, R., Lindsey, N. J., Biondi, B., & Ellsworth, W. L., 2021\.  Low-magnitude seismicity with a downhole distributed acoustic sensing array—examples from the forge geothermal experiment, Journal of Geophysical Research: Solid Earth, 126(1), e2020JB020462.
  * Li et al. (2018) Li, Z., Peng, Z., Hollis, D., Zhu, L., & McClellan, J., 2018.  High-resolution seismic event detection using local similarity for large-n arrays, Scientific reports, 8(1), 1646.
  * Li et al. (2021) Li, Z., Shen, Z., Yang, Y., Williams, E., Wang, X., & Zhan, Z., 2021.  Rapid response to the 2019 ridgecrest earthquake with distributed acoustic sensing, AGU Advances, 2(2), e2021AV000395.
  * Lindsey & Martin (2021) Lindsey, N. J. & Martin, E. R., 2021.  Fiber-optic seismology, Annual Review of Earth and Planetary Sciences, 49, 309–336.
  * Lindsey et al. (2019) Lindsey, N. J., Dawe, T. C., & Ajo-Franklin, J. B., 2019.  Illuminating seafloor faults and ocean dynamics with dark fiber distributed acoustic sensing, Science, 366(6469), 1103–1107.
  * Lindsey et al. (2020) Lindsey, N. J., Rademacher, H., & Ajo-Franklin, J. B., 2020.  On the broadband instrument response of fiber-optic das arrays, Journal of Geophysical Research: Solid Earth, 125(2), e2019JB018145.
  * LLC (2022) LLC, S., 2022.  Utah forge: Well 16a(78)-32 2022 stimulation silixa microseismic report, Tech. rep.
  * Maxwell et al. (2012) Maxwell, S., Raymer, D., Williams, M., & Primiero, P., 2012.  Tracking microseismic signals from the reservoir to surface, The Leading Edge, 31(11), 1300–1308.
  * McLennan (2022) McLennan, J., 2022.  Utah forge well 16a(78)-32 stimulation data (april, 2022).
  * Miah & Potter (2017) Miah, K. & Potter, D. K., 2017.  A review of hybrid fiber-optic distributed simultaneous vibration and temperature sensing technology and its geophysical applications, Sensors, 17(11), 2511.
  * Neidell & Taner (1971) Neidell, N. S. & Taner, M. T., 1971.  Semblance and other coherency measures for multichannel data, Geophysics, 36(3), 482–497.
  * Paitz et al. (2021) Paitz, P., Edme, P., Gräff, D., Walter, F., Doetsch, J., Chalari, A., Schmelzbach, C., & Fichtner, A., 2021.  Empirical investigations of the instrument response for distributed acoustic sensing (das) across 17 octaves, Bulletin of the Seismological Society of America, 111(1), 1–10.
  * Rossi et al. (2023) Rossi, C., Grigoli, F., Gasperini, P., Gandolfi, S., Cocorullo, C., Gukov, T., & Macini, P., 2023.  Estimation of amplitude noise reduction as a function of depth recorded by a deep vertical array (northern italy), Seismological Research Letter, p. to appear.
  * Shinohara et al. (2022) Shinohara, M., Yamada, T., Akuhara, T., Mochizuki, K., & Sakai, S., 2022.  Performance of seismic observation by distributed acoustic sensing technology using a seafloor cable off sanriku, japan, Frontiers in Marine Science, p. 466.
  * Stork et al. (2020) Stork, A. L., Baird, A. F., Horne, S. A., Naldrett, G., Lapins, S., Kendall, J.-M., Wookey, J., Verdon, J. P., Clarke, A., & Williams, A., 2020.  Application of machine learning to microseismic event detection in distributed acoustic sensing data, Geophysics, 85(5), KS149–KS160.
  * Tuinstra et al. (2023) Tuinstra, K., Grigoli, F., Lanza, F., Rinaldi, A. P., Fichtner, A., & Wiemer, S., 2023.  Locating clustered seismicity using distance geometry solvers: applications for sparse and single-borehole das networks, arXiv preprint arXiv:2309.16317.
  * Wang et al. (2021) Wang, X., Zhan, Z., Williams, E. F., Herráez, M. G., Martins, H. F., & Karrenbach, M., 2021.  Ground vibrations recorded by fiber-optic cables reveal traffic response to covid-19 lockdown measures in pasadena, california, Communications Earth & Environment, 2(1), 160.
  * Wang et al. (2020) Wang, Z., Lu, B., Ye, Q., & Cai, H., 2020.  Recent progress in distributed fiber acoustic sensing with $\phi$-otdr, Sensors, 20(22), 6594.
  * Webster et al. (2013) Webster, P., Wall, J., Perkins, C., & Molenaar, M., 2013.  Micro-seismic detection using distributed acoustic sensing, in 2013 SEG Annual Meeting, OnePetro.
  * Zeng et al. (2022) Zeng, X., Bao, F., Thurber, C. H., Lin, R., Wang, S., Song, Z., & Han, L., 2022\.  Turning a telecom fiber-optic cable into an ultradense seismic array for rapid postearthquake response in an urban area, Seismological Society of America, 93(2A), 853–865.
  * Zhidong et al. (2019) Zhidong, C., Shize, W., Wei, L., Fei, L., Chong, W., Liuyi, M., & Qing, L., 2019\.  Application of walkaway-vsp based on joint observation by das and geophones in the tarim basin, northwest china, in SEG International Exposition and Annual Meeting, OnePetro.
  * Zhu et al. (2023) Zhu, W., Biondi, E., Li, J., Yin, J., Ross, Z. E., & Zhan, Z., 2023.  Seismic arrival-time picking on distributed acoustic sensing data using semi-supervised learning, arXiv preprint arXiv:2302.08747.