# Lattice simulation of $SU(2)$ dark glueball with machine learning

Min-Huan Chu  INPAC, Key Laboratory for Particle Astrophysics and Cosmology
(MOE), Shanghai Key Laboratory for Particle Physics and Cosmology, School of
Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China
Yang Yuanqing Scientific Computering Center, Tsung-Dao Lee Institute, Shanghai
Jiao Tong University, Shanghai 200240, China  Jun-Hui Lai
<EMAIL_ADDRESS>INPAC, Key Laboratory for Particle Astrophysics and
Cosmology (MOE), Shanghai Key Laboratory for Particle Physics and Cosmology,
School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China  Wei Wang<EMAIL_ADDRESS>INPAC, Key Laboratory for
Particle Astrophysics and Cosmology (MOE), Shanghai Key Laboratory for
Particle Physics and Cosmology, School of Physics and Astronomy, Shanghai Jiao
Tong University, Shanghai 200240, China  Southern Center for Nuclear-Science
Theory (SCNT), Institute of Modern Physics, Chinese Academy of Sciences,
Huizhou 516000, Guangdong Province, China  Jialu Zhang  INPAC, Key Laboratory
for Particle Astrophysics and Cosmology (MOE), Shanghai Key Laboratory for
Particle Physics and Cosmology, School of Physics and Astronomy, Shanghai Jiao
Tong University, Shanghai 200240, China  Qianteng Zhu  INPAC, Key Laboratory
for Particle Astrophysics and Cosmology (MOE), Shanghai Key Laboratory for
Particle Physics and Cosmology, School of Physics and Astronomy, Shanghai Jiao
Tong University, Shanghai 200240, China

###### Abstract

We study the mass and scattering cross section of $SU(2)$ glueballs as dark
matter candidates using lattice simulations. We employ both naive and improved
$SU(2)$ gauge actions in $3+1$ dimensions with several $\beta$ values, and
adopt both the tranditional Monte Carlo method and the flow-based model based
on machine learning techniques to generate lattice configurations. The mass of
the scalar glueball with $J^{PC}=0^{++}$ and the NBS wave function are
calculated. Using the Runge-Kutta method, we extract the glueball interaction
potential and scattering cross section. From the observational constraints, we
obtain the lower bound of the mass of scalar glueball candidates as potential
components of dark matter.

##  I Introduction

Dark matter (DM) that gives rise to about a quarter of the universe’s mass
remains one of the most mysterious objects in particle physics. Primarily
recognized through its gravitational influence on galaxies and cosmic
structures, it might be still undetectable through electromagnetic means. The
weakly interacting massive particles miracle (WIMP), the suggestive
coincidence of the weak coupling DM with a proper density that would have
thermally been generated, motivated experimental tests through their direct
scattering, decays to visible cosmic rays, and DM productions at collider
experiments, but no positive result has been reported so far (for an
incomplete list on this topics, please see Refs. Bai:2010hh ; Fox:2011fx ;
PandaX-II:2016vec ; Slatyer:2016qyl ; LUX:2016ggv ; ATLAS:2017bfj ; PandaX-
II:2017hlx ; PICO:2017tgi ; XENON:2018voc ; ATLAS:2021kxv ; PandaX-4T:2021bab
; Foster:2021ngm ; PandaX:2023toi ; PandaX:2022ood ; XENON:2023cxc ).

Aside from WIMP, several other interesting DM models exist. An interesting
scenario suggests dark matter is likely to be glueballs, particles composed of
strongly interacting dark gluons bound together without fermions, stemming
from a confining dark $SU(N)$ gauge theory extrapolated from Quantum
Chromodynamics (QCD). These glueballs in the dark sector are hypothesized to
interact with standard model particle mainly through gravitational forces and
may contribute significantly to the dark matter content Carlson:1992fn ;
Boddy:2014yra ; Kribs:2016cew ; Dienes:2016vei ; Carenza:2022pjd ;
Carenza:2023eua .

Glueballs, as hypothesized entities within QCD, are bound states of gluons,
characterized by the absence of quark constituents. These entities emerge from
non-Abelian gauge theories, notably $SU(N)$ lattice gauge theories, where non-
perturbative effects support the formation of such composite particles
Soni:2016gzf ; Soni:2016yes ; Acharya:2017szw . This feature surprisingly
bridges the gap between the domain of strong interactions and cosmological
phenomena, suggesting a novel view where the dynamics of the strong force
support the mysterious nature of dark matter. The similarities underscore the
appealing hypothesis that glueballs may indeed constitute the dark matter
component in these scenarios Carenza:2022pjd . Theoretical studies and
simulations within this framework have been important in substantiating the
glueball dark matter hypothesis, offering a robust statement that integrates
the microcosmic interactions of particle physics with the large-scale
structures and dynamics observed in the universe.

In nonperturbative studies of glueball in QCD, lattice calculations play a
crucial role Morningstar:1999rf ; APE:1987ehd ; Chen:2005mg ; Gui:2012gx ;
Athenodorou:2020ani ; Sun:2017ipk ; Gui:2019dtm . This non-perturbative method
offers detailed insights into key properties of glueballs, including the
spectrum and decay properties. Consequently, it reveals their potential
significance across a spectrum of physical phenomena. Recent advances in
lattice calculations Yamanaka:2019yek ; Yamanaka:2019aeq  have been also
crucial in determining the scattering cross section of dark glueballs within
lattice $SU(2)$ Yang-Mills theory.

On the other side, machine learning algorithms are increasingly employed to
generate lattice configurations, enabling more efficient and accurate
investigations of glueball properties and behaviors. This computational
approach allows for handling complex data sets and extracting insights that
might be challenging to distinguish through traditional methods Kanwar:2020xzo
; Boyda:2020hsi ; Cranmer:2023xbe . One of the promising methods in recent
years is the flow-based model Albergo:2019eim ; Boyda:2020hsi ;
Albergo:2021bna ; Abbott:2022zhs  (for a review, please see Refs.
Albergo:2021vyo ; Cranmer:2023xbe ), which is inspired by the concepts and
techniques of machine learning, especially deep generative models. A flow-
based model is a type of neural network that can learn a complicated
probability distribution from data, and generate new samples from it, without
using Markov chain Monte Carlo methods. Moreover, flow-based model can also
learn the potential features and correlations of the data, and provide
insights into the structure and dynamics of the system. Flow-based models
impose an invertible transformation constructed by neural networks
transforming a simple distirbution to a complex one to perform the maximum
likelihood estimation. Kullback-Leibler divergence is used to quantify the
differences between the two distributions Kullback:1951zyt . To relieve the
mode-collapse problem in this approach, one can adopt a quenching algorithm
Chen:2022ytr .

In the remnant of this work, we will use the conventional Monte Carlo method
and machine learning approach to generate lattice configurations for the
$SU(2)$ gauge theory with zero fermion flavors ($N_{f}=0$) and several
coupling constants $\beta$. Subsequently, we compare the mass of the $0^{++}$
glueball obtained from these distinct simulations. Furthermore, we calculate
the inter potential for $0^{++}$ glueball and determine the interaction
coupling constants in effective Lagrangian. Finally, we obtain the relation
between the glueball mass and scattering cross section. By confronting
theoretical results with experimental data, we extract the lower bound on the
mass of scalar glueball candidates as potential components of dark matter.

##  II Methodology

In this section, we will exhibit the theoretical basis of our lattice
simulations. This comprises the application of the $SU(N)$ gauge theory action
within both the Monte Carlo method and machine learning frameworks. We
elaborate on the specific machine learning methodology employed, focusing on
the flow-based model. Additionally, this section also includes the strategy of
investigate the dark glueball, displaying how to determine its mass, coupling
constants of effective Lagrangian, and scattering cross section.

###  II.1 $SU(N)$ lattice gauge theory

The $SU(N)$ gauge field theory was established for describing interactions of
non-Abelian gauge boson, with its $SU(3)$ variant being crucial for modeling
strong interactions. Explicitly one can discretize lattice to simulate QCD in
continuum with variable gauge links $U_{\mu}(n)$ with path integral approach
Wilson:1974sk . Standard simulations easily incorporate plaquettes
$U_{\mu\nu}(n)$, defined as:

$\displaystyle
U_{\mu\nu}(n)=U_{\mu}(n)U_{\nu}(n+\hat{\mu})U^{\dagger}_{\mu}(n+\hat{\nu})U^{\dagger}_{\nu}(n),$
(1)

with $n$ indicating the lattice sites, and $\mu,\nu$ representing gauge link
directions. From the continuum action with variable $A_{\mu}(n)$, the
discretized action is constructed with variable $U_{\mu}(n)$ as:

$\displaystyle
S[U]=\frac{\beta}{N}\sum_{n,\mu<\nu}\mathrm{Re}\;\mathrm{tr}\left[1-U_{\mu\nu}(n)\right],$
(2)

where $N$ denotes the color degree of freedom in the $SU(N)$ gauge field, and
$\beta$ correlates to the continuum coupling constant as $\beta=2/g^{2}$. This
formulation is known as the naive action.

To minimize lattice spacing discretization effects, the Lüscher-Weisz
(improved) action is commonly used Luscher:1984xn , including both plaquettes
and rectangles:

$\displaystyle
R_{\mu\nu}(n)=U_{\mu}(n)U_{\mu}(n+\hat{\mu})U_{\nu}(n+2\hat{\mu})$
$\displaystyle
U^{\dagger}_{\mu}(n+\hat{\mu}+\hat{\nu})U^{\dagger}_{\mu}(n+\hat{\nu})U^{\dagger}_{\nu}(n).$
(3)

The tree-level correction is constructed as:

$\displaystyle S[U]=\frac{\beta}{2N}\mathrm{Re}$
$\displaystyle\Big{[}\sum_{n,\mu<\nu}c_{0}\mathrm{tr}\left(1-U_{\mu\nu}(n)\right)$
$\displaystyle+$
$\displaystyle\sum_{n,\mu,\nu}c_{1}\mathrm{tr}\left(1-R_{\mu\nu}(n)\right)\Big{]},$
(4)

where the coefficients for plaquettes and rectangles are $c_{0}=5/3$ and
$c_{1}=-1/12$, respectively.

Besides, it is acknowledged that simulations targeting glueballs require
extensive statistics, necessitating a judicious choice in lattice volume; it
cannot be excessively large. Conversely, given the substantial mass of
glueballs, the lattice spacing should not be overly small. An effective
solution to this challenge is employing an anisotropic lattice, characterized
by a larger length and finer lattice spacing in the temporal direction, which
is conducive for glueball simulations. The naive action in such a scenario is
given by:

$\displaystyle
S[U]=\frac{\beta}{N}\Big{[}\frac{1}{\xi}\sum_{n,i<j}\mathrm{Re}\;\mathrm{tr}\left(1-U_{ij}(n)\right)$
$\displaystyle+\xi\sum_{n,\mu}\mathrm{Re}\;\mathrm{tr}\left(1-U_{0\mu}(n)\right)\Big{]},$
(5)

where $\xi$ denotes the anisotropy factor, typically selected as an integer
like $2,3,4,5,\ldots$. Morningstar and Peardon pioneeringly applied
anisotropic lattices to glueball calculations Morningstar:1997ff .

Our investigation is an attempt of applying the self-interacting gauge fields
to dark matter study, and thereby we choose $N=2$. We utilize both isotropic
and anisotropic naive and improved actions for Monte Carlo simulations and
naive action with machine learning simulations, which will be discussed in
Sec.III.

###  II.2 Flow-based model

In a flow based model, an invertible and differentiable function maps a prior
distribution $r(z)$ of $z$ on some manifold $X$ to the model distribution
$q(x)$ of $x$ on the same manifold to approximate the target density $p(x)$.
The model distribution is determined via the change-of-variable formula:

$\displaystyle
q(x)=r(z)[J(z)^{-1}]=r(z)\left|\underset{kl}{\det}\frac{\partial
f_{k}(z)}{\partial z_{l}}\right|^{-1},$ (6)

where the function $f$, or flow, is composed of coupling layers in the neural
network, and is well-designed to make the Jacobian tractable such that the
model density can be efficiently computed.

The prior distribution can be chosen as any simple and calculable form. In
this calculation, a uniform distribution with respect to the Haar measure of
$SU(2)$ is sampled to construct a prior distribution.

To ensure gauge equivariance, a link variable $U_{\mu}(n)$ is shifted to a
$1\times 1$ Wilson loop $U_{\mu\nu}(n)$, and the flow must be a transformation
acting on an element of $SU(2)$ gauge group and is equivariant under matrix
conjugationAbbott:2022hkm . After a flow transformation, the loop variable
should be changed back to the link to form a new gauge configuration. The flow
makes a loop $U_{\mu\nu}(n)\to U^{\prime}_{\mu\nu}(n)$ and the new link is
updated as:

$\displaystyle
U^{\prime}_{\mu}(n)=U^{\prime}_{\mu\nu}(n)U^{\dagger}_{\mu\nu}(n)U_{\mu}(n).$
(7)

The eigenvalues of an $SU(2)$ matrix A can be simply written as
$\lambda_{1}=e^{i\theta}$ and $\lambda_{2}=e^{-i\theta}$. The diagonalization
of A can be conducted as $XUX^{-1}$, where $U={\rm
diag}(\lambda_{1},\lambda_{2})$. The invertible flow $f:G\to G$ should be a
kernel, meaning that $f\left(XUX^{-1}\right)=XF\left(U\right)X^{-1}$ for all
$U,X\in G$. The kernel is to transform untraced loops and should move the
density between conjugacy classes while preserving structure within those
classes.

To make the Jacobian of the flow efficient to calculate, we choose a subset of
the loops to be active and let them be the input of the flow. The
corresponding updated links will affect some other loops, which are denoted as
passive. The rest of the loops are denoted as frozen and will be the
coefficients of the transformation of the active loops. In this construction,
the Jacobian is ensured to be a triangular matrix.

The flow on $SU(N)$ is chosen to be the neural spline flow. First, the loops
are splited into three parts,
$x=\left[x_{\mathrm{active}},x_{\mathrm{passive}},x_{\mathrm{frozen}}\right]$.
Then, generate the parameters $\mathbf{\theta}=NN\left(x_{frozen}\right)$,
where $NN$ is an arbitrary neural network. After fixing the parameters, one
can update the loops
$x^{\prime}_{\mathrm{active}}=f_{\theta}\left(x_{\mathrm{active}}\right)$.
Finally, we get
$x^{\prime}=\left[x^{\prime}_{\mathrm{active}},x^{\prime}_{\mathrm{passive}},x_{\mathrm{frozen}}\right]$.
The form of the neural spline flow is proposed in Durkan:2019nsq , and a
detailed description of the flow-based model can be found in Ref.
Albergo:2021vyo .

To relieve the mode-collapse problem which is common in flow-based models, we
apply a quenching algorithm in our approach. The gauge coupling locates in the
same place as $\hbar$ in the partition function, which suggests that $g\to 0$
will be a classical limit. For small $g$, or large $\beta$ (which is
equivalent to the inverse temperature), one pays a large price for field
configurations that do not minimize the action, and the model favours some
special modes. In contrast, when $g\to\infty$, $\beta$ goes to zero, the Yang-
Mills action disappears completely, and all field configurations are not
suppressed and contribute to the path integral. Thereby we train the model
from small $\beta$, and gradually increase $\beta$ to the aimed value.

###  II.3 Dark glueballs and scattering

####  II.3.1 Glueball Spectrum

The $0^{++}$ glueball operator in lattice simulations, organized according to
parity, spin, and charge conjugation, is given as:

$\displaystyle\hat{O}(t,\vec{x})=\mathrm{Re}\;\;\mathrm{tr}\sum_{\mu<\nu}U_{\mu\nu}(t,\vec{x}),$
(8)

with $U_{\mu\nu}(t,\vec{x})$ symbolizing a plaquette. To account for the
vacuum state sharing the $0^{++}$ glueball state’s quantum numbers, its
expectation value is subtracted, reformulating the operator as
$\hat{O}^{\prime}(t,\vec{x})=\hat{O}(t,\vec{x})-\langle\hat{O}(t,\vec{x})\rangle$.
We can use two-point correlation functions to determine the mass of glueball,

$\displaystyle C_{2}(t)=\sum_{\vec{x},\vec{y}}\langle
0|\hat{O}(t,\vec{x})\hat{O}(0,\vec{y})|0\rangle.$ (9)

For the large timeslice, the two-point correlation function is reduced to

$\displaystyle C_{2}(t)\sim Ce^{-E_{0}t}.$ (10)

The effective energy $E_{0}$ represents the energy of the lowest glueball
(ground state) and when the momentum is zero, it is reduced to the mass
$m_{g}$.

In contrast to the naive action, the operator formulation for the $0^{++}$
glueball in the context of the improved action incorporates an additional
rectangular term, as defined in Eq.(3). The operator for the improved action
is thus articulated as:

$\displaystyle\hat{O}_{\mathrm{imp}}(t,\vec{x})=\mathrm{Re}\mathrm{tr}\sum_{\mu<\nu}\left[U_{\mu\nu}(t,\vec{x})+\frac{1}{20}R_{\mu\nu}(t,\vec{x})\right].$
(11)

For subsequent analyses, we employ this refined operator representation for
the $0^{++}$ glueball within the framework of configurations derived from the
improved action.

####  II.3.2 The two body scattering cross section

The cross section of dark glueballs scattering can be inferred from the
interaction glueball potential in a non-relativistic limit. In our study, we
employ the Bonn approximation to determine the cross section from the internal
glueball potential. The $0^{++}$ internal glueball potential is modeled as a
spherically symmetric Yukawa and Gaussian potential:

$\displaystyle V_{\mathrm{Yukawa}}(r)=Y\frac{e^{-m_{g}r}}{4\pi r},$ (12)
$\displaystyle
V_{\mathrm{Gaussian}}(r)=Z^{(1)}e^{-\frac{(m_{g}r)^{2}}{8}}+Z^{(2)}e^{-\frac{(m_{g}r)^{2}}{2}}.$
(13)

where $m_{g}$ represents the mass of the $0^{++}$ glueball.

In our investigation, we employ both Yukawa and Gaussian potentials, focusing
on the spherically symmetric cross section which necessitates the
consideration of the $k\to 0$ limit. Utilizing the Bonn approximation, the
differential cross section is derived as follows:

$\displaystyle\frac{d\sigma}{d\Omega}=\frac{4m^{2}_{g}}{K^{2}}\left|\int_{0}^{\infty}rV(r)\sin(Kr)dr\right|^{2},$
(14)

where $K=2k\sin\frac{\theta}{2}$ represents the momentum transfer with $k$ and
$\theta$ being the magnitude and scattering angle of the relative momentum
between two interacting particles, respectively. $V(r)$ denotes the potential.

Alternatively, given the scalar effective theory of the $0^{++}$ glueball, and
considering the $\phi^{3}$ and $\phi^{4}$ interaction, the effective
Lagrangian for glueballs can be expressed as:

$\displaystyle\mathcal{L}=\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{1}{2}m^{2}\phi^{2}-\frac{\lambda_{3}}{3!}\phi^{3}-\frac{\lambda_{4}}{4!}\phi^{4}.$
(15)

Figure 1: The Feynman diagram of $t$ channel for scalar field elastic
scattering.

Corresponding to the contributions in the potential, the $\phi^{4}$
interaction term gives a delta function, while it is difficult to be obtained
numerically. Therefore, only $\phi^{3}$ interaction is considered for this
Lagrangian. Perturbative theory yields the differential cross section for two-
body elastic scattering $\phi(p_{1})\phi(p_{2})\to\phi(p_{3})\phi(p_{4})$. The
Feynman diagram in Fig.1 provides the $i\mathcal{M}$ matrix element:

$\displaystyle i\mathcal{M}=(-i\lambda_{3})^{2}\frac{i}{k^{2}+m^{2}}\times
4=\frac{-4i\lambda_{3}^{2}}{k^{2}+m_{g}^{2}}.$ (16)

Setting $k=p_{3}-p_{1}=0$, we obtain the differential cross section:

$\displaystyle\frac{d\sigma}{d\Omega}=\frac{|\mathcal{M}|^{2}}{64\pi^{2}(2m_{g})^{2}}=\frac{\lambda_{3}^{4}}{64\pi^{2}m_{g}^{6}}.$
(17)

The $\phi^{3}$ interaction generates Yukawa potential in the coordinate space
in the non-relativistic limit. If $V_{\mathrm{Yukawa}}(r)$ is inserted to
Eq.(14), the differential cross section is expressed as:

$\displaystyle\frac{d\sigma}{d\Omega}\Bigr{\|}_{\mathrm{Yukawa}}=\frac{Y^{2}}{4\pi^{2}m_{g}^{2}}.$
(18)

Comparing Eq.(17) and Eq.(18), one can find that the value of $\lambda_{3}$ is
determined by the coefficient $Y$:

$\displaystyle\lambda_{3}=2m_{g}\sqrt{|Y|}.$ (19)

To probe the $s$-wave behavior resulting from a spherically symmetric
potential, the wave splitting method is applied. This method is grounded in
the radial Schrödinger equation, in which the wave function asymptotically
approaches the following form as $r\to\infty$:

$\displaystyle\Psi(r)\xrightarrow{r\to\infty}\frac{A_{l}}{r}\sin\left[kr+\delta_{l}(k)\right].$
(20)

By solving the radial Schrödinger equation, taking into account the known form
and coefficients of the potential, we can extract the wave function $\Psi(r)$
and determine the phase shift $\delta_{l}(k)$. This enables the calculation of
both the differential and total cross sections for the $s$-wave in the $k\to
0$ limit:

$\displaystyle\frac{d\sigma_{s}}{d\Omega}=\lim_{k\to
0}\frac{1}{k^{2}}\sin^{2}\left[\delta(k)\right],\;\;\sigma\sim\sigma_{s}=4\pi\frac{d\sigma_{s}}{d\Omega}.$
(21)

Accordingly, with the values of $Y$ or $Z^{(1)}$, $Z^{(2)}$ established, one
is well-positioned to determine both the $s$-wave differential cross section
$d\sigma_{s}/d\Omega$ and the total cross section $\sigma$.

####  II.3.3 Nambu-Bethe-Salpeter wave function

To determine the coefficient $Y$ in the Yukawa potential or $Z^{(1)}$,
$Z^{(2)}$ in the Gaussian potential, we utilize the Nambu-Bethe-Salpeter (NBS)
wave function from lattice calculations. The NBS wave function, originally
proposed for hadrons Aoki:2013cra ; Ishii:2012ssm , is naturally applicable to
glueballs. It is defined as:

$\displaystyle\Psi_{g}(\vec{r})=\sum_{\vec{x}}\left\langle
0\left|\hat{O}(0,\vec{x}+\vec{r})\hat{O}(0,\vec{x})\right|GG\right\rangle,$
(22)

where $\hat{O}$ is a glueball operator, and $|GG\rangle$ represents a hadron
state comprising two glueballs.

The NBS wave function related to the Schrödinger equation with the internal
$0^{++}$ glueball potential. The spatial Schrödinger equation with a non-local
potential is expressed as:

$\displaystyle\frac{1}{m_{g}}\nabla^{2}\Psi_{g}(\vec{r})=\int
d^{3}r^{\prime}U(\vec{r}^{\prime},\vec{r})\Psi_{g}(\vec{r}^{\prime}),$ (23)

where $U(\vec{r}^{\prime},\vec{r})$ approximates a local, central potential
$U(\vec{r}^{\prime},\vec{r})=V(r)\delta^{3}(\vec{r}^{\prime}-\vec{r})$. The
relation between $V(r)$ and $\Psi_{g}$ is:

$\displaystyle
V(r)=\frac{1}{m_{g}}\frac{\nabla^{2}\Psi_{g}(\vec{r})}{\Psi_{g}(\vec{r})}.$
(24)

Replacing $V(r)$ with the Eq.(12) or Eq.(13) allows the determination of the
coefficient $Y$ or $Z^{(1)}$, $Z^{(2)}$ from $\Psi_{g}(\vec{r})$. In lattice
simulations, correlation functions to extract the NBS wave function are
constructed as:

$\displaystyle C_{g}(t,\vec{r})=\sum_{\vec{x}}\left\langle
0\left|\hat{O}(t,\vec{x}+\vec{r})\hat{O}(t,\vec{x})\hat{S}(0,0)\right|0\right\rangle,$
(25)

with $\hat{S}$ as the glueball operator source. Then, the wave function is
given as

$\displaystyle\Psi_{g}(\vec{r})=\lim_{t\to\infty}\frac{C_{g}(t,\vec{r})}{C_{g}(t,\vec{r}=0)}.$
(26)

In principle a 2-body state could serve as the source $\hat{S}$. However,
considering the $0^{++}$ glueball, $|GG\rangle$ and $|G\rangle$ share the same
quantum numbers. Additionally, one cannot split 1-body source into spatially
separated correlators, unlike 2- and 3-body sources. Hence, we select
$\hat{S}=\hat{O}$.

##  III Numerical simulations and results

In this section, we present our numerical findings derived from lattice
simulations. This encompasses detailed settings for configuration generation,
the computed results for the $0^{++}$ glueball mass, along with the
determination of the coupling coefficient for the effective Lagrangian and the
glueball scattering cross section.

###  III.1 Lattice setup

In the configuration generation phase for proper lattice spacing, we have
varied the coupling constant $\beta$ within the set $\\{2.0,2.2,2.4,2.6\\}$
for both naive and anisotropic naive actions in our $SU(2)$ simulations, while
for the improved actions, we also use $\beta=3.7$.

Table 1: Lattice Setup for $SU(2)$ Monte Carlo Simulations: Initial point generation using random $SU(2)$ matrices, followed by a warm-up phase with heatbath method, and a production phase employing a combination of one heatbath cycle and multiple overrelaxation cycles. The asymmetry factor in the simulation, represented by $\xi$, is defined in Eq.(5). In the consideration of two point function fit, we choose $Nt=\xi\times Ns$ as volume. action | $\beta$ | $\xi$ | Volume | heat+over | Warm up | Conf
---|---|---|---|---|---|---
naive | 2.0 | 1 | $10^{4}$ | 1+2 | 500 | $1024\times 5$
naive | 2.0 | 3 | $8^{3}\times 24$ | 1+4 | 500 | $5000\times 5$
naive | 2.2 | 1 | $10^{4}$ | 1+3 | 700 | $1024\times 5$
naive(ML) | 2.2 | 1 | $10^{4}$ | / | / | $1024$
naive | 2.2 | 3 | $10^{3}\times 30$ | 1+4 | 700 | $3000\times 5$
naive | 2.4 | 2 | $16^{3}\times 32$ | 1+3 | 900 | $2000\times 4$
naive | 2.6 | 2 | $18^{3}\times 36$ | 1+5 | 1000 | $1000\times 6$
improved | 3.7 | 1 | $16^{4}$ | 1+5 | 1000 | $1500\times 6$

Table. 1 delineates the configurations employed in our study. During the
generation process, we have used both heatbath and overelaxation methods for
the $SU(2)$ gauge field. To enhance signal clarity, we implemented APE (Array
Processor with Emulator) smearing for the gauge links as described in
APE:1987ehd ,

$\displaystyle
U^{\mathrm{smear}}_{\mu}(n)=U_{\mu}(n)+\alpha\sum_{\pm\nu\neq\mu}U_{\nu}(n)U_{\mu}(n+\hat{\nu})U^{\dagger}_{\nu}(n+\hat{\mu}).$
(27)

In our simulations, we select a smearing parameter of $\alpha=0.1$. A total of
20 APE smearing iterations are uniformly applied across all configurations.
Employing the established relationship between the string tension $\sigma$ and
the scale parameter $\Lambda$ as documented in Allton:2008ty , we determine
the lattice spacing in terms of $\Lambda^{-1}$ for each simulation setting
Bali:2000gf . The computed lattice spacings are listed in Table 2.

In this study, for data stability, in the calculation of glueball mass, we use
basic isotropic configurations with $\beta$ values of $2.0$ and $2.2$. We also
includes a configuration based on machine learning at $\beta=2.2$, and
anisotropic configurations at $\beta=2.4$ and $2.6$. At the same time, for the
NBS wave function calculation, we use anisotropic basic configurations with a
wider range of $\beta$ values, such as $2.0$, $2.2$, $2.4$, and $2.6$, along
with a isotropic naive action configuration with $\beta=2.2$ from machine
learning. For all calculations, we also apply an improved isotropic action
configuration at $\beta=3.7$.

Table 2: The lattice spacings in string tension $\sigma$ unit and $SU(2)$ scale $\Lambda$ for configurations of several $\beta$ based on naive and improved action. action | $\beta$ | $a\sqrt{\sigma}$ | $a[\Lambda^{-1}]$
---|---|---|---
naive | 2.0 | 0.655(10) | 0.3837(58)
naive | 2.2 | 0.5045(57) | 0.2954(33)
naive | 2.4 | 0.3331(61) | 0.1950(36)
naive | 2.6 | 0.1718(19) | 0.1006(11)
improved | 3.7 | 0.2244(74) | 0.1314(43)

The preceding section details the configuration setup employed in the Monte
Carlo methodology. In contrast, for the machine learning based approach, we
restrict our simulations to the naive isotropic action at $\beta=2.2$ with a
lattice volume of $10^{4}$. Using configurations generated by the machine
learning approach, we subsequently apply a very few times of Monte Carlo
filtering processes to perform multiple screening iterations.

Figure 2: Comparison of computer costs from Monte Carlo and machine learning
methods.

In this investigation, we also undertake a comparative analysis of the
computational time consumption between the Monte Carlo method and a machine
learning approach. When implementing the Monte Carlo calculation for $10^{4}$
lattice configurations, the processing time is approximately 0.17 seconds per
configuration. This implies that the initial warm-up process requires around
170 seconds, and the total time for generation, including both warm-up and
production, is calculated as (Conf$\times$0.17+1700) seconds. In contrast, the
machine learning method involves distinct training and production phases,
consuming 21 hours for training and Conf$\times$0.12 seconds for production.
Fig.2 illustrates a cost comparison between the Monte Carlo and machine
learning methods. Based on this analysis, it can be inferred that the Monte
Carlo method is more efficient for small-scale computations, whereas machine
learning demonstrates greater suitability for extensive-scale simulations.
Therefore, configurations derived from machine learning exhibit reduced
correlation and demand less computational time. From this perspective, machine
learning emerges as an efficient method for generation.

However, the implementation of machine learning is not without challenges. One
significant issue is the substantial memory requirement during the training
process, which may exceed the capacity of limited-core systems. Additionally,
the persistent problem of mode-collapse in machine learning necessitates the
continued use of the Monte Carlo method for configuration filtering. Given
these considerations, our primary simulations employ the Monte Carlo approach,
while machine learning is utilized for a comparative analysis.

###  III.2 $0^{++}$ glueball mass

The effective mass of the $0^{++}$ glueball in the framework of $SU(2)$ gauge
theory is ascertainable through the analysis of the two-point correlation
function $C_{2}(t)$, as delineated in Eq.(9). Utilizing the glueball operator
definition provided in Eq.(11), we select $\mu$ and $\nu$ along the spatial
directions to characterize the $0^{++}$ glueball. The correlation data are
subsequently fitted to an exponential decay in Euclidean time. In light of the
periodic boundary conditions applied during the configuration generation, the
functional form of $C_{2}(t)$ is modeled as:

$\displaystyle C_{2}(t)\sim C_{0}e^{-m_{g}t}+C_{0}e^{-m_{g}(N_{t}-t)}.$ (28)

The fitting results are illustrated in Fig. 3. The reduced chi-squared
statistic, $\chi^{2}/\text{d.o.f.}$, suggests an effective fit with the dual
exponential terms, underscoring the reliability of our derived effective mass
estimates, which are enumerated in Table 3. Moreover, we extrapolate the
$0^{++}$ glueball mass to the continuum limitAllton:2008ty , and determine the
results $m_{g}=7.31(22)\Lambda$. In parallel, our simulations using the
improved action suggest an effective glueball mass of $m_{g}=7.78(37)\Lambda$.
When we compare these results with those obtained using the naive action, we
find that the improved action, even at such large lattice spacing
$a=0.1314(43)\Lambda$, has smaller discrete effect than naive action
configurations. We have put together all these results in Tab. 3, which also
shows very similar values for the naive action at $\beta=2.2$, achieved
through both machine learning and Monte Carlo methods. This comparison
confirms that conventional Monte Carlo methods and Machine Learning techniques
are consistent with each other. These results from different methods underline
the usefulness of Machine Learning in particle physics, especially for
complicated tasks like calculating the glueball mass.

Table 3: The results of effective masses of $0^{++}$ glueball are shown in lattice unit and in $\Lambda$. action | $\beta$ | $am_{g}$ | $m_{g}[\Lambda]$
---|---|---|---
naive | 2.0 | 1.685(93) | 4.39(25)
naive | 2.2 | 1.364(55) | 4.66(20)
naive(ML) | 2.2 | 1.331(55) | 4.50(19)
naive | 2.4 | 0.885(14) | 4.54(11)
naive | 2.6 | 0.763(16) | 7.58(18)
naive | continuum limit | / | 7.31(22)
improved | 3.7 | 1.022(34) | 7.78(37)

Figure 3: This plots show $C_{2}(t)$ results and fit situations, where the
$am_{g}$ and $chi^{2}/d.o.f.$ in the fit are explicitly given.

###  III.3 Scattering cross section

In the previous section, we address the determination of the glueball
interaction potential using the NBS wave function, as outlined in Eq.(24).
Tackling the numerical second partial derivative in this context is
challenging. Given the potential’s spherical symmetry, we approach this by
solving a radial partial differential equation:

$\displaystyle\frac{d^{2}R(r)}{dr^{2}}+\frac{2}{r}\frac{dR(r)}{dr}-m_{g}V(r)R(r)=0,$
(29)

Here, $R(r)$ is the radial part of $\Psi_{g}(\vec{r})$. The solution’s angular
component is expressed as spherical harmonics with quantum numbers $l,m$,
leading to $\Psi(\vec{r})=R(r)Y_{l,m}(\theta,\phi)$. Averaging
$\Psi_{g}(\vec{r})$ over a spherical shell simplifies $Y_{l,m}(\theta,\phi)$
to $1$, effectively reducing the NBS wave function $\Psi_{g}(\vec{r})$ to
$R(r)$ times a constant.

To proceed, we calculate the one-body correlation function as stated in
Eq.(25) and average it over spherical shells, approximating it as $R(r,t)$. We
also examine its dependence on time $t$. Fig. 4 shows this at $\beta=2.4$ in
naive action Monte Carlo simulations. We observe that $R(r,t)$ changes very
little with large $t$, so we use the averaged results from these areas. Our
findings for $R(r)$ from different lattice setups are presented in Fig. 5,
showing similar patterns.

Figure 4: $t$ dependence of radial NBS wave function $R(r)$.

Figure 5: This figure shows the radial part of the NBS wave function across
various lattice setups. We normalize all data at $r=0.3837\Lambda$, and show
the chosen $t$ range in the middle of each panel. As noted in the main text,
the configurations labeled ‘naive’ are based on anisotropic lattices with
finer lattice spacing in the $t$ direction. In contrast, the results from the
improved action and those obtained using machine learning are based on
isotropic lattices.

Considering Yukawa and Gaussian forms of potential energy, we set boundary
conditions and adjust coefficients $Y$ or $Z_{1}$, $Z_{2}$ to find solutions
for $R(r,Y)$ or $R(r,Z_{1},Z_{2})$ using the fourth-order Runge-Kutta method.
We then measure the differences between the $R(r,Y)$ or $R(r,Z_{1},Z_{2})$
curves and the data of $R(r)$, selecting the curve that most closely matches
our data to determine the coefficients.

Figure 6: The plot shows the results by solving the differential equation with
Yukawa and Gaussian potential.

Fig. 6 compares the solutions with Yukawa and Gaussian potentials to our
lattice data. This shows that solving the differential equation Eq.(29) with
the Runge-Kutta method aligns well with our simulations, indicating reliable
results. We thus obtained the values for $Y$ in the Yukawa potential, and
$Z_{1}$, $Z_{2}$ in the Gaussian potential. For clarity, we reformulate these
equations as follows:

$\displaystyle V_{\mathrm{Yukawa}}(r)=-20.53(0.82)\frac{e^{-m_{g}r}}{4\pi r},$
(30) $\displaystyle V_{\mathrm{Gaussian}}(r)=9.07(0.29)\Lambda
e^{-\frac{(m_{g}r)^{2}}{8}}$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;-16.5(1.2)\Lambda
e^{-\frac{(m_{g}r)^{2}}{2}}.$ (31)

Our current analysis shows that the Gaussian potential better matches the
original data, as we see in Fig. 6. We use a combination of bands from several
lattice sets for each potential, helping us find the final coefficients in
these potentials. Looking at the relationship between the Yukawa potential’s
coefficient $Y$ and the coupling coefficient $\lambda_{3}$ in the scalar
glueball Lagrangian (see Eq.(19)), we find
$\lambda_{3}=2m_{g}\sqrt{Y}=66.2(2.4)\Lambda$.

Because of these findings, we decide to use both potentials in our cross
section results. We treat the difference between these two types of potential
as a way to understand the systematic uncertainty in what we find.

Having established the values of $Y$, $Z_{1}$, and $Z_{2}$, we are now
equipped to formulate the Schrödinger equation for the radial wave function
$\phi(r)$, which is related to $r\Psi(r)$, in the presence of a non-zero
transfer momentum $k$. This equation is expressed as:

$\displaystyle\left(\frac{d^{2}}{dr^{2}}+k^{2}-m_{g}V(r)\right)\phi(r)=0$ (32)

The wave splitting method suggests that $\phi(r)$ approaches an asymptotic
form, $\phi(r)\propto\sin(kr+\delta(k))$, at large distances. Utilizing the
Runge-Kutta numerical method, we can accurately solve for $\phi(r)$. Applying
the parameters derived from the Yukawa and Gaussian potentials, we calculate
the $s$-wave total cross section as follows:

$\displaystyle\sigma=(1.13\sim 2.49)\Lambda^{-2}$ (33)

which contains the systematic and statistic uncertainties.

With the obtained results for $\sigma$ and $m_{g}$, we are now in a position
to elucidate the relationship between the ratio $\sigma/m_{g}$ and $m_{g}$.
Fig. 7 presents our findings alongside data from experimental
studiesRandall:2008ppe ; Peter:2012jh ; Harvey:2015hha ; Wittman:2017gxn . It
is evident from our analysis that the ratio $\sigma/m_{g}$ tends to decrease
as $m_{g}$ increases. Additionally, according to experimental data, it is
estimated that the value of $m_{g}$ is likely around $0.4$GeV.

Figure 7: This graph shows the ratio $\sigma/m_{g}$ against the glueball mass
$m_{g}$. It incorporates the constraints from several research groups that
have reported on the cross section of self-interacting dark matter. The band
in the graph represents the range encompassed by our calculations. For context
and comparison, the graph also includes results from established literature:
$\sigma/m=0.7\mathrm{cm}^{2}/\mathrm{g}$ as in Ref. Randall:2008ppe ,
$\sigma/m=0.1\mathrm{cm}^{2}/\mathrm{g}$ from Ref. Peter:2012jh ,
$\sigma/m=0.47\mathrm{cm}^{2}/\mathrm{g}$ from Ref. Harvey:2015hha , and
$\sigma/m=2.0\mathrm{cm}^{2}/\mathrm{g}$ from Ref. Wittman:2017gxn .

##  IV Summary and outlook

In this work, we have conducted an analysis of $0^{++}$ glueballs within the
$SU(2)$ gauge theory using both Monte Carlo simulations and machine learning
techniques. We have used these two approaches to generate various lattice
configurations, and accordingly extracted the glueball mass using these
configurations based on naive and improved actions. Then the glueball
interaction potential, crucial for extracting the interaction coupling
constant in effective quantum field theory and determining the dark glueball
scattering cross section, is extracted. We then established a connection
between the scattering cross section ($\sigma$) and the glueball mass
($m_{g}$), a vital aspect for understanding glueball behavior in various
physical scenarios. Our findings show that the ratio $\sigma/m_{g}$ decreases
with increasing $m_{g}$, aligning with values determined in experimental
studies. This correlation, along with estimated glueball mass $m_{g}\approx
0.3\;\mathrm{GeV}$ from experimental data, has significant implications for
dark matter research.

The application of machine learning alongside traditional Monte Carlo
simulations not only showcases the potential of computational advancements in
theoretical physics but also paves the way for more refined and efficient
research methodologies in the future. Our results can contribute to the
broader understanding of glueball dynamics and their interactions in the
universe Kang:2022jbg , particularly in relation to dark matter. This study
encourages further exploration into complex particle systems and their
interactions, which could be pivotal in understanding the mysteries of dark
matter and its fundamental forces.

## Acknowledgments

We are very grateful to Ying Chen, Longcheng Gui, Zhaofeng Kang, Hang Liu,
Jianglai Liu, Peng Sun, Jinxin Tan, Lingxiao Wang, Yibo Yang, Haibo Yu, Ning
Zhou, Jiang Zhu for valuable discussions. This work is supported in part by
Natural Science Foundation of China under grant No. 12125503 and 12335003. The
computations in this paper were run on the Siyuan-1 cluster supported by the
Center for High Performance Computing at Shanghai Jiao Tong University and
Advanced Computing East China Subcenter.

##  Appendix A Set Lattice Spacing from Wilson Loop

In the realm of lattice gauge theory, the expectation value of the Wilson loop
is crucial. We use this value to set the lattice spacing for each lattice
configuration. The static interaction potential for the $SU(2)$ gauge field is
determined as described in Bali:2000gf :

$\displaystyle
V(n_{r})=\lim_{t\to\infty}\frac{1}{a}\ln\frac{W(n_{r},n_{t})}{W(n_{r},n_{t}+1)},$
(34)

where $W(n_{r},n_{t})$ represents the Wilson loop’s expectation value with
length $n_{r}$ and width $n_{t}$ in lattice units. At short distances, this
potential behaves as $1/n_{r}$, and at larger distances, it exhibits a linear
growth, represented by the string tension $\sigma$.

One can express $aV(n_{r})$ as a function of $n_{r}$:

$\displaystyle aV(n_{r})=A+\frac{B}{n_{r}}+a^{2}\sigma n_{r}.$ (35)

For each lattice configuration, we fit the Wilson loop data to this formula to
obtain the string tension $a^{2}\sigma$. Additionally, using the relationship
between string tension and the energy scale $\Lambda$ for $SU(2)$ gauge theory
Teper:1998kw :

$\displaystyle\Lambda=0.586(45)\sqrt{\sigma},$ (36)

we finally determine the lattice spacing in units of $\Lambda$.

##  Appendix B Determining Wave Function from a Specified Potential

In Sec. III.3, we outline our method to determine coefficients for a specific
potential using the radial wave function $R(r)$, satisfying:

$\displaystyle\frac{d^{2}R(r)}{dr^{2}}+\frac{2}{r}\frac{dR(r)}{dr}-m_{g}V(r)R(r)=0.$
(37)

Here, $V(r)$ may represent Yukawa or Gaussian potentials, each requiring one
or two coefficients. Our approach is to use a test coefficient in $V(r)$ and
solve the differential equation via the Runge-Kutta method, starting from
$R(r_{\mathrm{ini}})$ and $R^{\prime}(r_{\mathrm{ini}})$. The initial
condition is derived from the NBS wave function of lattice data, replacing
differentiation with difference.

Applying the fourth-order Runge-Kutta method, we rewrite the above equation
as:

$\displaystyle
R^{\prime\prime}(r)=f(r,R,R^{\prime})=\frac{2}{r}R^{\prime}(r)-m_{g}R(r)V(r).$
(38)

Then we can obtain the solution $R(r)$ by iterating the following equations
with step length $h$:

$\displaystyle k_{1}R=hr,\;\;k_{1}R^{\prime}=hf(r,R,R^{\prime}),$ (39)
$\displaystyle k_{2}R=h(R^{\prime}+\frac{1}{2}k_{1}R^{\prime}),$ (40)
$\displaystyle
k_{2}R^{\prime}=hf(r+\frac{1}{2}h,R+\frac{1}{2}k_{1}R,R^{\prime}+\frac{1}{2}k_{1}R^{\prime}),$
(41) $\displaystyle k_{3}R=h(R^{\prime}+\frac{1}{2}k_{2}R^{\prime}),$ (42)
$\displaystyle
k_{3}R^{\prime}=hf(r+\frac{1}{2}h,R+\frac{1}{2}k_{2}R,R^{\prime}+\frac{1}{2}k_{2}R^{\prime}),$
(43) $\displaystyle k_{4}R=h(R^{\prime}+k_{3}R^{\prime}),$ (44) $\displaystyle
k_{4}R^{\prime}=hf(r+h,R+k_{3}R,R^{\prime}+k_{3}R^{\prime}),$ (45)
$\displaystyle r_{1}=r+h,$ (46) $\displaystyle
R_{1}=R+\frac{1}{6}(k_{1}R+2k_{2}R+2k_{3}R+k_{4}R),$ (47) $\displaystyle
R^{\prime}_{1}=R^{\prime}+\frac{1}{6}(k_{1}R^{\prime}+2k_{2}R^{\prime}+2k_{3}R^{\prime}+k_{4}R^{\prime}),$
(48)

where $r_{1}$, $R_{1}$, and $R^{\prime}_{1}$ are the iteratively updated
values at $r+h$. Using test coefficients in $V(r)$, we approximate the
corresponding solution, denoted as $R_{\mathrm{test}}(r)$. We then assess the
distance between this test solution and our lattice data. Then by varying the
coefficients within a certain range and finding the minimal distance, we
eventually identify the corresponding coefficients for $V(r)$.

## References

  * (1) Y. Bai, P. J. Fox and R. Harnik, JHEP 12, 048 (2010) doi:10.1007/JHEP12(2010)048 [arXiv:1005.3797 [hep-ph]].
  * (2) P. J. Fox, R. Harnik, J. Kopp and Y. Tsai, Phys. Rev. D 84, 014028 (2011) doi:10.1103/PhysRevD.84.014028 [arXiv:1103.0240 [hep-ph]].
  * (3) A. Tan et al. [PandaX-II], Phys. Rev. Lett. 117, no.12, 121303 (2016) doi:10.1103/PhysRevLett.117.121303 [arXiv:1607.07400 [hep-ex]].
  * (4) T. R. Slatyer and C. L. Wu, Phys. Rev. D 95, no.2, 023010 (2017) doi:10.1103/PhysRevD.95.023010 [arXiv:1610.06933 [astro-ph.CO]].
  * (5) D. S. Akerib et al. [LUX], Phys. Rev. Lett. 118, no.2, 021303 (2017) doi:10.1103/PhysRevLett.118.021303 [arXiv:1608.07648 [astro-ph.CO]].
  * (6) M. Aaboud et al. [ATLAS], JHEP 01, 126 (2018) doi:10.1007/JHEP01(2018)126 [arXiv:1711.03301 [hep-ex]].
  * (7) X. Cui et al. [PandaX-II], Phys. Rev. Lett. 119, no.18, 181302 (2017) doi:10.1103/PhysRevLett.119.181302 [arXiv:1708.06917 [astro-ph.CO]].
  * (8) C. Amole et al. [PICO], Phys. Rev. Lett. 118, no.25, 251301 (2017) doi:10.1103/PhysRevLett.118.251301 [arXiv:1702.07666 [astro-ph.CO]].
  * (9) E. Aprile et al. [XENON], Phys. Rev. Lett. 121, no.11, 111302 (2018) doi:10.1103/PhysRevLett.121.111302 [arXiv:1805.12562 [astro-ph.CO]].
  * (10) G. Aad et al. [ATLAS], Phys. Rev. D 103, no.11, 112006 (2021) doi:10.1103/PhysRevD.103.112006 [arXiv:2102.10874 [hep-ex]].
  * (11) Y. Meng et al. [PandaX-4T], Phys. Rev. Lett. 127, no.26, 261802 (2021) doi:10.1103/PhysRevLett.127.261802 [arXiv:2107.13438 [hep-ex]].
  * (12) J. W. Foster, M. Kongsore, C. Dessert, Y. Park, N. L. Rodd, K. Cranmer and B. R. Safdi, Phys. Rev. Lett. 127, no.5, 051101 (2021) doi:10.1103/PhysRevLett.127.051101 [arXiv:2102.02207 [astro-ph.CO]].
  * (13) X. Ning et al. [PandaX], Nature 618, no.7963, 47-50 (2023) doi:10.1038/s41586-023-05982-0
  * (14) D. Zhang et al. [PandaX], Phys. Rev. Lett. 129, no.16, 161804 (2022) doi:10.1103/PhysRevLett.129.161804 [arXiv:2206.02339 [hep-ex]].
  * (15) E. Aprile et al. [XENON], Phys. Rev. Lett. 131, no.4, 041003 (2023) doi:10.1103/PhysRevLett.131.041003 [arXiv:2303.14729 [hep-ex]].
  * (16) E. D. Carlson, M. E. Machacek and L. J. Hall, Astrophys. J. 398, 43-52 (1992) doi:10.1086/171833
  * (17) K. K. Boddy, J. L. Feng, M. Kaplinghat and T. M. P. Tait, Phys. Rev. D 89, no.11, 115017 (2014) doi:10.1103/PhysRevD.89.115017 [arXiv:1402.3629 [hep-ph]].
  * (18) G. D. Kribs and E. T. Neil, Int. J. Mod. Phys. A 31, no.22, 1643004 (2016) doi:10.1142/S0217751X16430041 [arXiv:1604.04627 [hep-ph]].
  * (19) K. R. Dienes, F. Huang, S. Su and B. Thomas, Phys. Rev. D 95, no.4, 043526 (2017) doi:10.1103/PhysRevD.95.043526 [arXiv:1610.04112 [hep-ph]].
  * (20) P. Carenza, R. Pasechnik, G. Salinas and Z. W. Wang, Phys. Rev. Lett. 129, no.26, 26 (2022) doi:10.1103/PhysRevLett.129.261302 [arXiv:2207.13716 [hep-ph]].
  * (21) P. Carenza, T. Ferreira, R. Pasechnik and Z. W. Wang, Phys. Rev. D 108, no.12, 12 (2023) doi:10.1103/PhysRevD.108.123027 [arXiv:2306.09510 [hep-ph]].
  * (22) A. Soni and Y. Zhang, Phys. Rev. D 93, no.11, 115025 (2016) doi:10.1103/PhysRevD.93.115025 [arXiv:1602.00714 [hep-ph]].
  * (23) A. Soni and Y. Zhang, Phys. Lett. B 771, 379-384 (2017) doi:10.1016/j.physletb.2017.05.077 [arXiv:1610.06931 [hep-ph]].
  * (24) B. S. Acharya, M. Fairbairn and E. Hardy, JHEP 07, 100 (2017) doi:10.1007/JHEP07(2017)100 [arXiv:1704.01804 [hep-ph]].
  * (25) C. J. Morningstar and M. J. Peardon, Phys. Rev. D 60, 034509 (1999) doi:10.1103/PhysRevD.60.034509 [arXiv:hep-lat/9901004 [hep-lat]].
  * (26) M. Albanese et al. [APE], Phys. Lett. B 192, 163-169 (1987) doi:10.1016/0370-2693(87)91160-9
  * (27) Y. Chen, A. Alexandru, S. J. Dong, T. Draper, I. Horvath, F. X. Lee, K. F. Liu, N. Mathur, C. Morningstar and M. Peardon, et al. Phys. Rev. D 73, 014516 (2006) doi:10.1103/PhysRevD.73.014516 [arXiv:hep-lat/0510074 [hep-lat]].
  * (28) L. C. Gui et al. [CLQCD], Phys. Rev. Lett. 110, no.2, 021601 (2013) doi:10.1103/PhysRevLett.110.021601 [arXiv:1206.0125 [hep-lat]].
  * (29) A. Athenodorou and M. Teper, JHEP 11, 172 (2020) doi:10.1007/JHEP11(2020)172 [arXiv:2007.06422 [hep-lat]].
  * (30) W. Sun, L. C. Gui, Y. Chen, M. Gong, C. Liu, Y. B. Liu, Z. Liu, J. P. Ma and J. B. Zhang, Chin. Phys. C 42, no.9, 093103 (2018) doi:10.1088/1674-1137/42/9/093103 [arXiv:1702.08174 [hep-lat]].
  * (31) L. C. Gui, J. M. Dong, Y. Chen and Y. B. Yang, Phys. Rev. D 100, no.5, 054511 (2019) doi:10.1103/PhysRevD.100.054511 [arXiv:1906.03666 [hep-lat]].
  * (32) N. Yamanaka, H. Iida, A. Nakamura and M. Wakayama, Phys. Rev. D 102, no.5, 054507 (2020) doi:10.1103/PhysRevD.102.054507 [arXiv:1910.07756 [hep-lat]].
  * (33) N. Yamanaka, H. Iida, A. Nakamura and M. Wakayama, Phys. Lett. B 813, 136056 (2021) doi:10.1016/j.physletb.2020.136056 [arXiv:1910.01440 [hep-ph]].
  * (34) G. Kanwar, M. S. Albergo, D. Boyda, K. Cranmer, D. C. Hackett, S. Racanière, D. J. Rezende and P. E. Shanahan, Phys. Rev. Lett. 125, no.12, 121601 (2020) doi:10.1103/PhysRevLett.125.121601 [arXiv:2003.06413 [hep-lat]].
  * (35) D. Boyda, G. Kanwar, S. Racanière, D. J. Rezende, M. S. Albergo, K. Cranmer, D. C. Hackett and P. E. Shanahan, Phys. Rev. D 103, no.7, 074504 (2021) doi:10.1103/PhysRevD.103.074504 [arXiv:2008.05456 [hep-lat]].
  * (36) K. Cranmer, G. Kanwar, S. Racanière, D. J. Rezende and P. E. Shanahan, Nature Rev. Phys. 5, no.9, 526-535 (2023) doi:10.1038/s42254-023-00616-w [arXiv:2309.01156 [hep-lat]].
  * (37) M. S. Albergo, G. Kanwar and P. E. Shanahan, Phys. Rev. D 100, no.3, 034515 (2019) doi:10.1103/PhysRevD.100.034515 [arXiv:1904.12072 [hep-lat]].
  * (38) M. S. Albergo, G. Kanwar, S. Racanière, D. J. Rezende, J. M. Urban, D. Boyda, K. Cranmer, D. C. Hackett and P. E. Shanahan, Phys. Rev. D 104, no.11, 114507 (2021) doi:10.1103/PhysRevD.104.114507 [arXiv:2106.05934 [hep-lat]].
  * (39) R. Abbott, M. S. Albergo, D. Boyda, K. Cranmer, D. C. Hackett, G. Kanwar, S. Racanière, D. J. Rezende, F. Romero-López and P. E. Shanahan, et al. Phys. Rev. D 106, no.7, 074506 (2022) doi:10.1103/PhysRevD.106.074506 [arXiv:2207.08945 [hep-lat]].
  * (40) M. S. Albergo, D. Boyda, D. C. Hackett, G. Kanwar, K. Cranmer, S. Racanière, D. J. Rezende and P. E. Shanahan, [arXiv:2101.08176 [hep-lat]].
  * (41) S. Kullback and R. A. Leibler, The Annals of Mathematical Statistics 22, no.1, 79-86 (1951) doi:10.1214/aoms/1177729694
  * (42) S. Chen, O. Savchuk, S. Zheng, B. Chen, H. Stoecker, L. Wang and K. Zhou, Phys. Rev. D 107, no.5, 056001 (2023) doi:10.1103/PhysRevD.107.056001 [arXiv:2211.03470 [hep-lat]].
  * (43) K. G. Wilson, Phys. Rev. D 10, 2445-2459 (1974) doi:10.1103/PhysRevD.10.2445
  * (44) M. Luscher and P. Weisz, Commun. Math. Phys. 98, no.3, 433 (1985) [erratum: Commun. Math. Phys. 98, 433 (1985)] doi:10.1007/BF01205792
  * (45) C. J. Morningstar and M. J. Peardon, Phys. Rev. D 56, 4043-4061 (1997) doi:10.1103/PhysRevD.56.4043 [arXiv:hep-lat/9704011 [hep-lat]].
  * (46) R. Abbott, M. S. Albergo, A. Botev, D. Boyda, K. Cranmer, D. C. Hackett, G. Kanwar, A. G. D. G. Matthews, S. Racanière and A. Razavi, et al. PoS LATTICE2022, 036 (2023) doi:10.22323/1.430.0036 [arXiv:2208.03832 [hep-lat]].
  * (47) C. Durkan, A. Bekasov, I. Murray and G. Papamakarios, [arXiv:1906.04032 [stat.ML]].
  * (48) S. Aoki, N. Ishii, T. Doi, Y. Ikeda and T. Inoue, Phys. Rev. D 88, no.1, 014036 (2013) doi:10.1103/PhysRevD.88.014036 [arXiv:1303.2210 [hep-lat]].
  * (49) N. Ishii et al. [HAL QCD], Phys. Lett. B 712, 437-441 (2012) doi:10.1016/j.physletb.2012.04.076 [arXiv:1203.3642 [hep-lat]].
  * (50) C. Allton, M. Teper and A. Trivini, JHEP 07, 021 (2008) doi:10.1088/1126-6708/2008/07/021 [arXiv:0803.1092 [hep-lat]].
  * (51) G. S. Bali, Phys. Rept. 343, 1-136 (2001) doi:10.1016/S0370-1573(00)00079-X [arXiv:hep-ph/0001312 [hep-ph]].
  * (52) S. W. Randall, M. Markevitch, D. Clowe, A. H. Gonzalez and M. Bradac, Astrophys. J. 679, 1173-1180 (2008) doi:10.1086/587859 [arXiv:0704.0261 [astro-ph]].
  * (53) A. H. G. Peter, M. Rocha, J. S. Bullock and M. Kaplinghat, Mon. Not. Roy. Astron. Soc. 430, 105 (2013) doi:10.1093/mnras/sts535 [arXiv:1208.3026 [astro-ph.CO]].
  * (54) D. Harvey, R. Massey, T. Kitching, A. Taylor and E. Tittley, Science 347, 1462-1465 (2015) doi:10.1126/science.1261381 [arXiv:1503.07675 [astro-ph.CO]].
  * (55) D. Wittman, N. Golovich and W. A. Dawson, Astrophys. J. 869, no.2, 104 (2018) doi:10.3847/1538-4357/aaee77 [arXiv:1701.05877 [astro-ph.CO]].
  * (56) Z. Kang, J. Zhu and J. Guo, Phys. Rev. D 107, no.7, 076005 (2023) doi:10.1103/PhysRevD.107.076005 [arXiv:2211.09442 [hep-ph]].
  * (57) M. J. Teper, [arXiv:hep-th/9812187 [hep-th]].


# Using 3D Shadows to Detect Object Hiding Attacks on Autonomous Vehicle
Perception

Zhongyuan Hau  Department of Computing
Imperial College London
<EMAIL_ADDRESS>Soteris Demetriou  Department of Computing
Imperial College London
<EMAIL_ADDRESS>Emil C. Lupu  Department of Computing
Imperial College London
<EMAIL_ADDRESS>
###### Abstract

Autonomous Vehicles (AVs) are mostly reliant on LiDAR sensors which enable
spatial perception of their surroundings and help make driving decisions.
Recent works demonstrated attacks that aim to hide objects from AV perception,
which can result in severe consequences. 3D shadows, are regions void of
measurements in 3D point clouds which arise from occlusions of objects in a
scene. 3D shadows were proposed as a physical invariant valuable for detecting
spoofed or fake objects. In this work, we leverage 3D shadows to locate
obstacles that are hidden from object detectors. We achieve this by searching
for void regions and locating the obstacles that cause these shadows. Our
proposed methodology can be used to detect an object that has been hidden by
an adversary as these objects, while hidden from 3D object detectors, still
induce shadow artifacts in 3D point clouds, which we use for obstacle
detection. We show that using 3D shadows for obstacle detection can achieve
high accuracy in matching shadows to their object and provide precise
prediction of an obstacle’s distance from the ego-vehicle.

###### Index Terms:

Autonomous Vehicles, LiDAR, 3D Object Detection, Obstacle Detection,
Automotive Security

##  I Introduction

Autonomous vehicles (AVs) are increasingly being deployed on public roads.
AV’s perception of the environment they operate in is fundamental to their
autonomy, guiding their driving decisions for safe and reliable operation.
These vehicles are commonly equipped with LiDAR sensors, which collect high
definition depth measurements stored in 3D point clouds. Objects detected from
the 3D point clouds are used by AVs to map the obstacles in the driving
environment, which is vital to the safety of the autonomous vehicle, its
passengers and surroundings. Object hiding attacks have been explored in
various works [1, 2, 3, 4] which demonstrate that it is possible to generate
adversarial objects that can evade object detection. Object detection is a
safety-critical function and failing to detect an object can lead to fatal
collisions.

Object hiding attacks. There are currently 2 classes of Object Hiding Attacks.
The first type is when the adversarial operations are performed on the
physical (target) object itself. One proposed method [1] is to place
adversarial objects on top of a target vehicle that evade point-cloud based
object detectors, achieving attack success rate of 80%. The adversarial
objects can be generated in a white-box or a black-box method. For the white-
box attack, the adversarial object is generated using a gradient-based
approach to minimize the confidence score of the target object (vehicle);
whereas the black-box attack chooses adversarial objects using a genetic
algorithm approach to iterate and improve on adversarial object meshes. Cao et
al. in [2] proposed using an optimization-based approach to generate
adversarial objects that evade detection by 3D object detectors. In [3], the
authors exploited the weakness of deep neural-networks (DNN) based detectors
and proposed an optimization-based method to generate physically realizable
adversarial objects that target multi-sensor fusion, effectively evading both
Camera and LiDAR based detectors.

The second type of object hiding attacks is when adversarial operations are
performed on the sensing modality itself. In [4], the authors proposed using
the LiDAR spoofing approach (Object Removal Attack, ORA-Random) to inject
random points in the target’s object bounding box to displace points in the
original object’s point cloud, causing mis-detection of the target object.
They showed that ORA was effective in damaging the performance of object
detectors for front-near objects and works well for both large objects such as
cars and small objects such as pedestrians and cyclists.

Motivation. Whilst many works demonstrate different ways to hide objects from
detection by object detectors, there is a lack of research into how we can
reliably detect such object hiding attacks. Cao et al. in [3] showed that
model-level defenses (i.e. input transformations and adversarial training) on
multi-sensor fusion models to increase their adversarial robustness does not
effectively decrease the attack success rates. The authors have also suggested
fusing more sensor modalities to increase the robustness, but argue that this
approach does not fundamentally address the weakness of using DNN-based object
detectors and would still be exploitable. To this end, we develop an object
hiding attack detection methodology that leverages a strong physical invariant
of a single sensing modality (LiDAR) to provide an orthogonal and reliable
detection of objects in a 3D scene, which is robust against object hiding
attacks.

Figure 1: Object Hiding Attacks from: (a) Tu et al.[1] (b) Cao et al.[5] (c)
Cao et al.[3] (d) Hau et al.[4]. Shadows of adversarial objects circled in
orange.

Our work. 3D shadows in LiDAR point clouds, introduced in [6], are a physical
phenomenon that is caused by occlusion of LiDAR laser pulses by objects in a
scene. The authors leverage these shadow artifacts as a physical invariant to
detect LiDAR spoofing attacks. We observed in Figure 1, that under object
hiding attacks, the objects hidden by the adversary in the scene, are not
detected by a DNN-based object detector, but still exist and occlude LiDAR
pulses to leave behind shadow artifacts in the 3D point cloud. We leverage
this observation and propose an obstacle detection methodology that searches a
region of interest in a 3D point cloud for objects based on the shadow
artifacts found in the 3D point cloud scene. Our proposed methodology of using
shadows provides an alternate approach to detecting objects in a 3D scene and
can be used as an orthogonal defense against object hiding attacks.

##  II Threat Model

In this work, we investigate the detection of unidentified objects using their
shadow artifacts in a 3D point cloud. Locating unidentified objects can help
us detect Object Hiding Attacks, a class of attacks on AV perception that aims
to cause mis-detection of objects in the sensed environment.

Threat model. We assume an adversary $\mathcal{A}$, that has the ability to
perform any of the state-of-the-art object hiding attacks. For Class 1 Object
Hiding Attacks, there are no perturbations to the 3D LiDAR point cloud, and
thus, the observation that object shadows exist holds. For Class 2, Object
Removal Attack (ORA-Random) [4] generate adversarial objects by perturbing the
point cloud. However, the point perturbations are limited to a short distance
of point shifting along the ray direction and do not affect the object’s
shadow.

Emulating object hiding attacks. Object Hiding Attacks are thus observed to
not alter the 3D shadows of the adversarial objects. Additionally, we also
assume an adversary that does not perturb the shadows. As such, to emulate
such object hiding attack for the evaluation of our proposed detection
mechanism, we simply ignore that these objects are detected by an object
detector and remove the labelled object from any output. In short, we assume
that the attacker has successfully conducted an object hiding attack and the
output labels do not contain the target object.

##  III Using 3D Shadows to Detect Hidden Objects

In this work, we use the key observation that 3D shadows are caused by
occlusion of LiDAR laser pulses by opaque objects. Adversarial objects in a 3D
point cloud scene that are subjected to Object Hiding Attacks, would not be
detected by an DNN-based object detector, but due to the physics of LiDAR, the
phenomenon of 3D shadows would still exist. The existence of the objects’ 3D
shadows in a 3D point cloud can be exploited to locate and identify
adversarial objects that evade object detection.

3D Shadows as physical invariant. 3D shadows have been introduced by Hau et
al. in [6] as a phyiscal invariant to verify genuine objects in a 3D point
cloud scene and detect object spoofing attacks. In object spoofing attacks [7,
8, 5, 9], point cloud measurements are injected into a scene to spoof objects
that are then erroneously detected by 3D-object detectors used by AVs. Shadow-
Catcher [4] detects such attacks by using the bounding box of objects
generated by object detectors to identify a shadow region behind the detected
object and inspects that shadow region for presence of point cloud
measurements, which are indicative of a spoofing attack. The authors have
shown that 3D shadows are a strong physical invariant that can be used to
verify genuine objects. In this work, we explore searching the scene for void
regions and attributing them to genuine objects.

3D Shadows to detect object hiding attacks. From this observation that 3D
shadows are strong a physical invariant that can be used to verify genuine
objects in a 3D scene, we propose a methodology that identifies void regions
in a 3D scene, attributes it to objects and corroborates this with the output
of object detectors to: 1) verify genuine objects and 2) detect potential
object hiding attacks in a scene.

##  IV Obstacle Detection Design

We propose a methodology that searches a region of interest in the 3D point
cloud scene for void regions, and analyzes them to either verify that the void
region is caused by the occlusion of a detected object or an unidentified
obstacle that could be an adversarial hidden object. This approach is
different from Shadow-Catcher [6] as Shadow-Catcher uses the bounding box of a
detected object to obtain the shadow region for analysis to detect spoofed
objects. With the Object Hiding Attack adversary we consider that the attacker
aims to hide the object and thus there is no bounding box for Shadow-Catcher
to perform analysis on – Shadow-Catcher is blind to an object hiding
adversary. As our approach is unable to identify the objects (i.e. classify
and label the object), in the rest of the paper, we call the objects detected
by our approach as obstacles. Objects hidden by an adversary will still have
shadows in a 3D point cloud, and thus can still be detected as an obstacle by
our methodology.

Identifying shadow regions. We first define a region of interest (RoI) to
analyse for shadows. We look into the front-near region of the ego-vehicle and
define a region of width 10m and length 30m in front of it, where we analyse
for shadows and subsequently identify obstacles. Next, from the RoI, a thin
volume of the point cloud around the ground level is extracted for analysis.
The RoI volume is discretized into cubic cells with width of 0.3m, and the
occupancy of these cells by LiDAR point measurements is checked. Cells that
are empty are clustered using DBSCAN [10] into shadow clusters.

The parameters used such as RoI size and cubic cell width can be tuned based
on the visibility requirements (i.e. a larger RoI for greater visibility) and
the resolution of the LiDAR (a LiDAR that uses less lasers would have a lower
resolution and thus, requires a larger cell width for analysis).

Attributing shadows to their obstacles and detected object. In a simple case
where objects in a scene are sparsely located, each shadow cluster would
correspond to a single obstacle. However, in crowded scenes, the shadows of
the objects overlap resulting in a single large cluster. In order to
distinguish the different obstacles that contribute to the shadow, we use the
center of the cells in each shadow cluster to generate a collection of
frustums to the LiDAR on the ego-vehicle (origin coordinates of the point
cloud). The points in the collection of frustums are then checked to determine
whether they fall within bounding boxes of objects identified by 3D Object
Detectors. The remaining points that do not fall in the bounding box of
detected objects are then further processed to identify and localize the
obstacles.

Localizing unidentified obstacles. The remaining points (that do not belong to
any identified objects) from the frustums of shadows are further processed as
follows. The points undergo clustering with DBSCAN where we get clusters of
points that according to the physics of LiDAR, contribute to the formation of
shadows and are indicative of an opaque obstacle in the scene. Finally, we
generate individual bounding boxes for each cluster that enclose their points.

###  IV-A Example 1: Single Object in Front-Near Region

We first demonstrate the use of the proposed obstacle detection in a simple
case where the scene (Fig. 2(a)) has a single Car object in the front-near
region of ego-vehicle, and the birds-eye-view (BEV) of the 3D point cloud in
Fig. 2(b).

(a) Image of Scene with 1 Object (Car) In Front-Near Region

(b) BEV of 3D point cloud with 1 object in Front-Near Region

Figure 2: Scene with 1 Object (Car) In Front-Near Region

We use our proposed Obstacle Detection methodology and first search for void
regions in the scene by discretizing the RoI into cells and checking if the
cells are occupied by points. The empty cells are clustered and the result is
shown in Fig. 3. Frustums are then generated (Fig. 4(a) from the LiDAR to the
cells in the shadow clusters, where the points in the frustums are clustered
as obstacles (Fig. 4(b)). Finally, the clusters of points, identified as
obstacles in the scene, are enclosed with bounding boxes (Fig. 4(c)).

Figure 3: Clustering of empty cells into shadow clusters – 1 Shadow Cluster.
The origin (0,0) is the LiDAR unit on the ego-vehicle.

(a) Frustums from Shadow Clusters to LiDAR

(b) Clusters of points in frustums

(c) Bounding Boxes of Obstacles identified from clusters

Figure 4: Processing of Shadow Clusters into Obstacles

Our proposed detection methodology was able to accurately detect the car in-
front of the ego-vehicle as an obstacle and in addition, picked up two other
obstacles which are the road divider and the traffic light post.

###  IV-B Example 2: Multiple Objects in Front-Near Region

Next, we demonstrate the effectiveness of the proposed obstacle detection
methodology in a scene where there are multiple objects in the front-near
region of the ego-vehicle. The scene is shown in Fig. 5, where there are
multiple cars parked along the side of the road.

Figure 5: Image of Scene with Multiple Objects In Front-Near Region.

Using the Obstacle Detection methodology, the scene is analysed for void
regions in the front-near area of interest and the void cells are clustered
together to form shadow clusters. Frustums from the LiDAR to the shadow
clusters are then generated (Fig. 6(a), the points in the frustums are
clustered to obtain point clusters of obstacles (Fig. 6(b) and their bounding
boxes are also generated (Fig. 6(c)).

(a) Frustums from Shadow Clusters to LiDAR

(b) Clusters of points in frustums

(c) Bounding Boxes of Obstacles identified from clusters

Figure 6: Processing of Shadow Clusters into Obstacles

From Fig. 6(c), the obstacles identified using the proposed Obstacle Detection
Methodology match well with the objects that are found in the scene (i.e. the
cars that are parked along the sides of the road).

###  IV-C Summary

We have proposed an Obstacle Detection Methodology that first searches an area
of interest in the front-near region of the ego-vehicle for void spaces that
are indicative of shadow regions. Ray frustums from the LiDAR to the shadow
clusters are generated and points that falls within these frustums are further
clustered to obtain point clusters of obstacles that are found in the scene.
From the demonstration of the 2 examples, we show that our proposed
methodology is capable of accurately localizing obstacles in a scene. Next, we
perform experiments to quantitatively measure the accuracy and effectiveness
of the proposed methodology in detecting obstacles.

##  V Experiments & Evaluation

Models & Datasets. We perform our evaluation of the effectiveness and
efficiency of our shadow detection methodology and object verification on the
KITTI dataset [11], which contains LiDAR measurements from scenes captured in
the real-world. Object Hiding Attacks are emulated by removing the target
object’s label from KITTI ground truth labels.

Evaluation scenarios. In our evaluation of the effectiveness and utility of
the proposed object hiding attack detection methodology, we design experiments
that aim to answer the following research questions:

RQ1: How well does the shadow identification and verification work in a benign
scenario? In this experiment, we use benign scenes from KITTI and evaluate how
well the shadow identification and verification procedure by measuring the
accuracy of matching shadows to ground truth labelled objects. For each
identified shadow, we sample the discretized space of the shadow and obtain
the frustum from the LiDAR to the shadow cell. Next, the frustum is analysed
for point measurements, if there are points in the frustum, the points are
corroborated with ground truth labels to determine if they belong to a
labelled object. We measure the accuracy (True Positives) of the shadow
identification and object verification procedure by counting the total matched
objects. We define False Negatives as labelled objects that we have missed,
and False Positives as detected obstacles that do not match the labelled
objects.

RQ2: How well does the methodology identify regions in a scene where an object
is subjected to Object Hiding Attack? We simulate Object Hiding Attacks and
evaluate how well we can localize the regions in a scene where objects exists
but are undetected by 3D Object Detectors. We use the procedure to identify
shadow regions and clusters of shadows. For each shadow cluster, their
discretized shadow space is sampled and the frustum from the LiDAR is
generated to analyse for points in these frustums. If the points in the
frustums do not belong to any detected objects, the un-labelled points are
recorded. Then, all the un-labelled points are clustered to form clusters of
un-labelled objects. We generate bounding boxes that enclose these clusters of
points and effectively use them to localize the object. We then match the un-
labelled objects to adversarially removed objects (from the KITTI ground
truth) and calculate the Birds-Eye-View 2D Intersection-over-Union (IoU) of
the ground truth object bounding boxes with those generated by localizing
objects from their shadows. We measure the average IoU of all the hidden
objects in the front-near region of interest from the ego-vehicle to evaluate
the accuracy of the localization approach. Additionally, we measure the
maximum distance of the nearest edge (to the ego-vehicle) of the bounding
boxes.

###  V-A Shadow Identification and Object Attribution

We run the proposed obstacle detection on 7480 scenes of the KITTI dataset.
For every shadow detected, the frustums are generated and points in frustums
are corroborated with ground truth bounding boxes of objects by checking if
these points fall in the labelled bounding boxes. In the 7480 scenes, there
are 10426 objects that are in the front-near region of interest. Out of these,
10259 (TPR 98.4%) were matched, showing a high accuracy in detecting objects.
We analyzed the missed objects (FNR 1.6%), and we found that the objects
missed by the obstacle detection are due to these objects being located too
close to the periphery of the RoI such the object is in the ROI but their
shadows fall outside the analysis region. The current width of the RoI is
chosen to be about 1.5$\times$ the width of a vehicle lane from the center of
the ego-vehicle and hence the objects missed are at least 1 vehicle lane away,
as the objects are closer to the center of the RoI, the shadows will be picked
up and the obstacle would be detected. The width of the RoI is a parameter
that can be increased if more visibility is required. We also recorded a False
Positive Rate of 11.9%, these are obstacles that are detected by the proposed
approach but do not match the ground-truth labelled objects. Upon further
analysis of these cases, we attribute the FPs to objects that are found
commonly in the environment such as road dividers and traffic lights, as
demonstrated in Section IV-A. Although the current approach is unable to
identify obstacles and provide object labels, it is able to faithfully detect
obstacles in a driving scene.

###  V-B Localizing Object Hiding Attack Detection

Next, we investigate how well the proposed obstacle detection is able to
localize objects and compare the localization to ground truth labelled
bounding boxes. In the 7480 scenes, obstacles were localized using the
proposed obstacle detection and we measure the 2D BEV IoU of the predicted
bounding box with the ground truth bounding boxes of objects in the front-near
region. The average IoU for the bounding boxes of obstacles detected in the
scenes is 33.2%. The low average IoU can be attributed to the predicted
bounding boxes only enclosing the point clusters found in the frustums (Fig.
7). The smaller predicted bounding boxes are due to "missing points" that can
arise from the DBSCAN clustering that label points that supposedly belonging
to the object but are too far away as noise, and hence are not included into
the clusters.

(a) BEV of Bounding Boxes (Predicted and Ground Truth)

(b) Another view of Bounding Boxes (Predicted and Ground Truth)

Figure 7: IoU example with a scene with multiple objects.

However, we observed that low IoU does not indicate that the predicted
bounding boxes are not well localized, rather they are just smaller in size.
As such, we used another metric where we measure the distance between the
nearest edges (to the ego-vehicle) of the predicted and ground truth bounding
boxes. This provides a metric of how well the obstacle detection predicts the
distance of the obstacle from the ego-vehicle. We found that the average
distance between the predicted and ground truth bounding box edges to be 1.8m
(standard deviation = 2.4m), which shows that the prediction of the obstacle
is very close to ground truth labels in terms of distance from the ego-
vehicle.

###  V-C Runtime Efficiency

We measure the runtime of the proposed obstacle detection using a machine
equipped with an Intel Core i7 Six Core Processor i7-7800X (3.5GHz) and 32GB
RAM. The average runtime for the full obstacle detection on 7480 scenes is
36.5 s/scene with a standard deviation of 7.7 s/scene. The majority of the
processing time for the obstacle detection pipeline is spent on the two DBSCAN
operations, first for clustering empty cells and second for clustering points
in the frustums. As such, the current implementation is not ideal for use in
real-time obstacle detection for autonomous driving decisions.

###  V-D Summary

We have evaluated the utility and efficiency of the proposed obstacle
detection methodology. We were able to detect 3D shadows and match them to the
obstacles that caused them with a high accuracy (98.4%). Although the proposed
obstacle detection is unable to predict the correct size of the object, it can
be used to estimate the obstacle’s distance away from the ego-vehicle. Our
current prototype implementation of our proposed obstacle detection mechanism
was found to be not ideal for real-time processing of AV scenes and we are
planning to work on optimisations to improving its efficiency.

##  VI Discussion

Detection approach. Our proposed obstacle detection approach uses shadows as
physical invariant of obstacles in a driving scene to locate these obstacles
in a scene. An approach used by Autoware [12] is to perform euclidean
clustering on the points in a down-sampled point-cloud to detect objects,
where its object detection result can be verified against the 3D object
detector’s output and a rule-based criteria can be used to determine if an
object is being hidden. While augmenting Autoware’s approach for detection is
straight-forward, our approach of using shadows as an orthogonal defense is
still useful as it provides a new perspective (looking at a physical
invariant) and opportunities for better interpreting anomalies.

Inability to identify type of object. Compared to DNN-based 3D object
detection models, our proposed obstacle detection is unable to recognize the
type of the object (i.e. unable to label the object). This is due to the use
of clustering algorithms to cluster points that are in close proximity to each
other to form clusters that represent the objects. The use of clustering
approach does not have any notion of how various object types are represented
in point cloud and hence, unable to predict the type of object. The clustering
approach that group points by proximity also leads to the shortcoming of poor
prediction of an object’s size. Future work could look into how different
modalities or machine learning approaches can be combined to enrich the
information from shadows and detected obstacles for improving the semantics of
our predictions.

Real-time processing requirements. Processing driving scenes in real-time is a
critical requirement for autonomous driving decision making. Due to the
computational overhead of our proposed approach, it is not suitable for real-
time obstacle detection. However, there is still merit for the current
approach to be used as an offline analysis tool that can aid in the detection
of anomalous objects. Future work can be done to improve on the efficiency of
the approach by exploring other clustering algorithms and localization
strategies to optimise the performance. Another possible approach would be to
have a lightweight implementation that identify an obstacle’s location by only
using the shadow clusters, doing away with the need of clustering points in
frustums to generate the bounding boxes.

Shadow removal adversary. The threat model we consider in this work assumes
that the adversary only performs object hiding attacks and does not take steps
to conceal the attack against our proposed obstacle detection. An adaptive
adversary can try to evade obstacle detection by performing LiDAR spoofing
attacks to inject LiDAR point measurements into the void region of "hidden
objects" to perturb their shadows. In future work, we plan to investigate the
robustness of using shadows to detect obstacles against a stronger adversary
that can simultaneously hide objects and perturb their shadows.

##  VII Conclusion & Future Work

In this paper, we leveraged the key observation that objects from object
hiding attacks leave behind shadow artifacts in their LiDAR 3D point cloud and
this can be exploited to detect objects hidden from 3D object detection. To
this end, we propose an obstacle detection methodology that first searches the
front-near region of the ego-vehicle for void regions in the 3D point cloud
and localizing obstacles that causes these shadows. We evaluated the
performance of the obstacle detection methodology and found that it can match
shadows to object with an accuracy of 98.4%. It is also able to predict the
distance of the obstacles to the ego-vehicle with an average precision of
1.8m. In future work, we hope to improve further on the utility of obstacle
detection, looking into ways to increase the precision of the localization of
obstacles and improving the efficiency to meet real-time processing
requirements. We also plan to evaluate the robustness of using shadows for
obstacle detection against an adaptive attacker that aims to remove shadows to
evade detection.

## Acknowledgment

We would like to thank Qi Alfred Chen and the anonymous reviewers for their
valuable feedback on our work.

## References

  * [1] J. Tu, M. Ren, S. Manivasagam, M. Liang, B. Yang, R. Du, F. Cheng, and R. Urtasun, “Physically realizable adversarial examples for lidar object detection,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 13 716–13 725.
  * [2] Y. Cao, C. Xiao, D. Yang, J. Fang, R. Yang, M. Liu, and B. Li, “Adversarial objects against lidar-based autonomous driving systems,” _arXiv preprint arXiv:1907.05418_ , 2019.
  * [3] Y. Cao, N. Wang, C. Xiao, D. Yang, J. Fang, R. Yang, Q. A. Chen, M. Liu, and B. Li, “Invisible for both camera and lidar: Security of multi-sensor fusion based perception in autonomous driving under physical-world attacks,” in _2021 IEEE Symposium on Security and Privacy (SP)_. IEEE, 2021, pp. 176–194.
  * [4] Z. Hau, K. T. Co, S. Demetriou, and E. C. Lupu, “Object removal attacks on lidar-based 3d object detectors,” _arXiv preprint arXiv:2102.03722_ , 2021\.
  * [5] Y. Cao, C. Xiao, B. Cyr, Y. Zhou, W. Park, S. Rampazzi, Q. A. Chen, K. Fu, and Z. M. Mao, “Adversarial sensor attack on lidar-based perception in autonomous driving,” in _Proceedings of the 2019 ACM SIGSAC Conference on Computer and Communications Security_ , 2019, pp. 2267–2281.
  * [6] Z. Hau, S. Demetriou, L. Muñoz González, and E. C. Lupu, “Shadow-catcher: Looking into shadows to detect ghost objects in autonomous vehicle 3d sensing,” in _26th European Symposium on Research in Computer Security (ESORICS)_ , 2021.
  * [7] J. Petit, B. Stottelaar, M. Feiri, and F. Kargl, “Remote attacks on automated vehicles sensors: Experiments on camera and lidar,” _Black Hat Europe_ , vol. 11, p. 2015, 2015.
  * [8] H. Shin, D. Kim, Y. Kwon, and Y. Kim, “Illusion and dazzle: Adversarial optical channel exploits against lidars for automotive applications,” in _International Conference on Cryptographic Hardware and Embedded Systems_. Springer, 2017, pp. 445–467.
  * [9] J. Sun, Y. Cao, Q. A. Chen, and Z. M. Mao, “Towards robust lidar-based perception in autonomous driving: General black-box adversarial sensor attack and countermeasures,” in _29th USENIX Security Symposium (USENIX Security 20)_. USENIX Association, Aug. 2020, pp. 877–894. [Online]. Available: https://www.usenix.org/conference/usenixsecurity20/presentation/sun
  * [10] M. Ester, H.-P. Kriegel, J. Sander, X. Xu _et al._ , “A density-based algorithm for discovering clusters in large spatial databases with noise.” in _KDD_ , vol. 96, no. 34, 1996, pp. 226–231.
  * [11] A. Geiger, P. Lenz, C. Stiller, and R. Urtasun, “Vision meets robotics: The kitti dataset,” _International Journal of Robotics Research (IJRR)_ , 2013\.
  * [12] S. Kato, S. Tokunaga, Y. Maruyama, S. Maeda, M. Hirabayashi, Y. Kitsukawa, A. Monrroy, T. Ando, Y. Fujii, and T. Azumi, “Autoware on board: Enabling autonomous vehicles with embedded systems,” in _2018 ACM/IEEE 9th International Conference on Cyber-Physical Systems (ICCPS)_. IEEE, 2018, pp. 287–296.