# E-RAFT: Dense Optical Flow from Event Cameras

Mathias Gehrig Mario Millhäusler∗ Daniel Gehrig Davide Scaramuzza
Dept. Informatics, Univ. of Zurich and
Dept. of Neuroinformatics, Univ. of Zurich and ETH Zurich
Equal contribution

###### Abstract

We propose to incorporate feature correlation and sequential processing into
dense optical flow estimation from event cameras. Modern frame-based optical
flow methods heavily rely on matching costs computed from feature correlation.
In contrast, there exists no optical flow method for event cameras that
explicitly computes matching costs. Instead, learning-based approaches using
events usually resort to the U-Net architecture to estimate optical flow
sparsely. Our key finding is that the introduction of correlation features
significantly improves results compared to previous methods that solely rely
on convolution layers. Compared to the state-of-the-art, our proposed approach
computes dense optical flow and reduces the end-point error by 23% on MVSEC.
Furthermore, we show that all existing optical flow methods developed so far
for event cameras have been evaluated on datasets with very small displacement
fields with maximum flow magnitude of 10 pixels. Based on this observation, we
introduce a new real-world dataset that exhibits displacement fields with
magnitudes up to 210 pixels and 3 times higher camera resolution. Our proposed
approach reduces the end-point error on this dataset by 66%.

## Multimedia Material

Code: https://github.com/uzh-rpg/E-RAFT
Dataset: https://dsec.ifi.uzh.ch
Video: https://youtu.be/dN8fl7-XfNw

##  1 Introduction

Figure 1: Optical flow predictions from our proposed approach trained on our
novel optical flow dataset for event-based vision. The second and third row
visualize input event data and images of the same scene respectively.

Optical flow estimation is a fundamental computer vision task that informs
about movement in image space. It serves as a fundamental building block for
visual odometry [28] for robotics and AR/VR, autonomous driving [16], action
recognition [31], high dynamic range (HDR) imaging [18], and computational
videography [25]. Traditionally, this problem is formulated as finding pixel
correspondences between two images.

While earlier top-performing methods were mostly based on energy minimization
techniques [13], current state of the art is dominated by neural network-based
approaches [34].

Decades of research on image-based optical flow estimation have resulted in
impressive results on public benchmarks [23]. However, image-based methods
struggle in presence of motion blur [6] or over-saturated image regions due to
limited dynamic range of image sensors. A promising sensor for optical flow
estimation is the event camera.

Event cameras drastically differ from image-based sensors in their way of
acquiring data. In contrast to frame-based cameras, which capture images at
regular intervals, event cameras report per-pixel brightness changes as a
stream of asynchronous events. For a detailed survey on event cameras, we
refer to [9]. The distinct advantages of event cameras are high temporal
resolution (microseconds), high dynamic range ($>120$ dB), and no motion blur.
All three aspects are essential for high-quality optical flow estimation in
real-world applications. However, these advantages are tied to the
fundamentally different data format resulting from event cameras.

Effectively extracting motion from an asynchronous event stream is a non-
trivial task. Unlike images, which can be directly used as input to
convolutional neural networks, event data is sparse, irregular, and
asynchronous in nature. Prior work on optical flow estimation from events
proposed spatio-temporal plane-fitting [3], variational optimization based on
discretized spatio-temporal volumes [2] and, more recently, convolutional
architectures on preprocessed event data [40, 41, 38]. Similarly to image-
based optical flow, the current state of the art is based on data-driven
neural network approaches. These data-driven methods have in common that they
are only evaluated for _sparse_ optical flow prediction while relying on a
U-Net architecture [29] to extract motion features from events converted into
a tensor representation.

In contrast to event-based optical flow methods, modern image-based optical
flow approaches make use of cost volumes to extract discriminant features [35,
32, 36, 34]. Typically, a neural network computes image features at a lower
resolution, which are then used to compute feature correlations between the
source and the target descriptors. The feature correlations are summarized in
a cost volume that is processed with CNN layers [32, 36]. Teed et al. [34]
chose instead to extract features from the cost volume via an iterative update
scheme.

In contrast to prior work, we focus on _dense_ optical flow estimation for
event cameras and utilize cost volumes instead of using the U-Net
architecture. First, we show that a state-of-the-art method designed for
images can be adopted to work exceptionally well on event-data, outperforming
all previously proposed methods consistently. Then, we exploit the fact that
event data is sequential and propose a differentiable warm-starting technique
that initializes the flow estimate for the next time step. Differently from
previous work [34], we perform backprogagation through time via the
differentiable warm-starting module. Our experiments suggest that
incorporating this warm-starting technique in the training pipeline improves
the performance on sequential event data further. More importantly, our
proposed approach can be directly used for fine-tuning our non-recurrent model
for sequential data because it does not introduce any additional parameters.

Until now, data driven optical flow methods on events have solely been
evaluated on the MVSEC dataset [40]. Our analysis suggests that the difficulty
of this dataset is very limited. For example, 95% of ground truth flow on the
outdoor evaluation sequence [40] has a magnitude of less than 3 pixels (see
figure 3). Based on this observation, we extend a recently proposed real-world
driving dataset [11] with optical flow ground truth. This dataset features a
maximum flow magnitude of up to 210 pixels and an event camera resolution that
is three times higher than MVSEC. We train and evaluate a state-of-the-art
event-based method on this dataset and show that our proposed approach
outperforms existing approaches that have so far only been tested on MVSEC by
23% EPE (0.27 vs 0.35) and by 66% EPE (0.79 vs 2.32) on our proposed dataset.

In summary, our contributions are

  * •

The first optical flow approach for event cameras that is designed and
evaluated for dense optical flow estimation.

  * •

The first optical flow method for event cameras, which utilizes cost volumes
and introduces recurrency to incorporate temporal priors. Previous event-based
optical flow methods rely on plane fitting, variational optimization or pure
CNN layers to estimate optical flow.

  * •

A novel real-world optical flow dataset for event-based vision with
significantly larger pixel displacements, higher resolution than current
datasets, and higher ground truth quality.

  * •

Controlled experiments and evaluation on both MVSEC and our proposed dataset
that suggest that both cost volumes and temporal recurrency are powerful
concepts for event-based optical flow estimation outperforming purely
convolutional architectures.

##  2 Related Work

###  2.1 Learning-based Optical Flow Estimation on Images

Data-driven approaches for optical flow estimation has been pioneered by
FlowNet [8]. In the following years, optical flow estimation has witnessed
great progress through principled neural network design [39, 14, 15, 21, 32,
33, 36, 34]. As a result a few core ingredients for high performance optical
flow estimation have emerged: _(i)_ coarse-to-fine processing, _(ii)_ the use
of feature correlation volumes, and _(iii)_ feature warping [39, 14, 15, 21,
32, 33, 36]. Teed et al. [34] builds on these ideas and extends them by
constructing correlation volumes at multiple scales and introducing an
iterative approach to extract optical flow from them, reaching state-of-the-
art results on the Sintel[6] and KITTI[12] benchmarks.

While these methods work well in standard conditions, their performance is
intricately tied to the quality of the images used. These algorithms are thus
not able to cope with image degradations, that occurr in challenging
scenarios, such as in HDR scenes and during high-speed motion.

In this work, we seek to address this intrinsic limitation of current
approaches, by proposing the use of event cameras [20] which do not suffer
from motion blur or limited dynamic range.

###  2.2 Event-based Optical Flow

Algorithmic methods for event-based optical flow estimation can be summarized
in three categories. First, iterative asynchronous methods inspired by [4, 10]
by the well-known Lucas-Kanade algorithm [22]. The latter also leverages image
data for photometric feature tracking. Second, plane fitting methods [5, 24]
that exploit the local plane-like shape of spatio-temporal event streams.
Third, variational optimization based approaches [2, 26] that incorporate
image data [26] or simultaneously estimate image intensity [2]. These methods
directly work on event data and often incorporate strong assumptions, such as
brightness constancy [2, 26], into their models. In contrast, our approach
learns abstractions and priors from data and does not rely on such
assumptions.

Another line of research pursues learning-based approaches and exclusively
adopts the self-supervised learning framework [41, 19, 37, 27, 40].
Supervision is either provided by images [40, 41, 19] or events only [37, 27].
However, these methods, except for [37], choose to mask their optical flow
prediction for evaluation due to limiting accuracy in regions where no events
are present. On the other hand, Ye et al. [37] assume that the scene is static
and predict optical flow via depth and ego-motion. All aforementioned
approaches have in common that they use the U-Net architecture to directly
estimate and evaluate _sparse_ optical flow from event representations. In
contrast to above methods, we predict and evaluate _dense_ optical flow, make
explicit use of longer temporal history, and leverage correlation features
from cost volumes to estimate pixel correspondences for large displacements.
Finally, we depart from the common theme of using a U-Net architecture for
event-based optical flow estimation while at the same time outperforming
existing methods consistently on multiple datasets.

###  2.3 Datasets for Event-based Optical Flow

Dataset | Train/Test Split | Groundtruth Quality | Resolution (MP) | Translational Motion | Outdoor
---|---|---|---|---|---
DVSFLOW16 | ✘ | high | 0.05 | ✘ | ✘
DVSMOTION20 | ✘ | high | 0.1 | ✘ | ✘
MVSEC | (✘) | medium/low | 0.1 | ✔ | ✔
Ours | ✔ | high | 0.3 | ✔ | ✔
Table 1: Comparison of existing event camera datasets for optical flow
estimation.

To the best of our knowledge, there are currently only three event camera
datasets that contain optical flow ground truth data. Table 1 summarizes and
relates them to our novel proposed dataset.

DVSFLOW16 [30] and DVSMOTION20 [1] are two small datasets designed for
evaluating event-based optical flow algorithms for rotational camera motion in
indoor scenes. Due to their small size and restricted camera motion, they are
not well suited for training and evaluation of neural network based
approaches.

More closely related to our work is the MVSEC dataset [40]. It contains 5
driving sequences and 4 indoor sequences in the same environment. Groundtruth
optical flow data is computed from lidar odometry and mapping for the outdoor
scenes. However, there are three major limitations of this dataset. First, the
ground truth quality is impaired by the lack of proper handling of occlusions
and moving objects. Second, the resolution of the event cameras is low with
only 0.1 Megapixels leading to small pixel displacements. We analyze this
aspect in detail in section 4. Finally, MVSEC does not introduce a clear train
and test split [40]. For example, the authors evaluate their proposed method
on the outdoor day 1 sequence but also train on this sequence to evaluate on
the indoor sequences. We address all these shortcomings by proposing a novel
dataset free of ambiguities from occlusions and artifacts from moving objects,
3 times higher camera resolution, and a well-defined training and test split
with established metrics to foster reproducibility.

##  3 Approach

Our approach takes inspiration from recent work on two-frame optical flow
estimation [34]. However, instead of images, we use consecutive packets of
events, which we convert to tensor-like event representations that are
compatible with standard CNNs. Two CNN encoders with shared weights create a
feature embedding to compute an all-pairs correlation volume at $1/8$-th of
the original resolution. The same encoder architecture, but with a different
set of weights, is used to encode context features for further processing.
Both outputs are used in an iteration scheme to extract per pixel flow, as
proposed in [34]. We then extend this formulation to a temporally recurrent
architecture, that processes sequences of event packets, by leveraging
differentiable warmstarting of optical flow (Sec. 3.3).

###  3.1 Problem Definition

Dense optical flow estimation is commonly defined as the task of finding pixel
correspondences between two images $I_{i}$ and $I_{j}$ [23, 6] denoted as flow
$F_{i\to j}$. For event cameras, we are interested in the analog setting of
finding pixel correspondence $F_{t_{i}\to t_{j}}$ between two timestamps
$t_{i}$ and $t_{j}$ given events in between both timestamps. In other words,
our setting is strictly sequential in contrast to image-based optical flow
estimation. For the sake of concise notation, the remaining text refers to
$F_{t_{i}\to t_{j}}$ as $F_{i\to j}$.

Different from frame-based cameras that acquire images at a specific rate,
event cameras respond to brightness changes in the scene asynchronously and
independently for each pixel. They do this by maintaining a per-pixel
reference log intensity. As soon as a sufficient deviation from this memorized
value is detected, an event is triggered and the memorized value is updated.
An event $e_{k}(t)=(x_{k},y_{k},t,p_{k})$ is a tuple containing information
about its pixel location $(x_{k},y_{k})$, the exact trigger time $t$ as well
as the 1-bit polarity $p_{k}$ that indicates whether the brightness change was
positive or negative. We define a sequence of events between time $t_{i}$ and
$t_{j}$ as ${\cal E}(t_{i},t_{j}):=\langle e(t)|t\in[t_{i},t_{j}]\rangle$.

###  3.2 RAFT for Event Data

Previous approaches for event-based optical flow estimation with neural
networks build an input representation given event data between $t_{i}$ and
$t_{i+1}=t_{i}+\Delta t$ and directly regresses optical flow $F_{i\to i+1}$
[40, 41, 27]. Instead, we propose to regress $F_{i\to i+1}$ from _two_ short
consecutive event sequences ${\cal E}(t_{i-1},t_{i})$ and ${\cal
E}(t_{i},t_{i+1})$.

Our baseline approach builds on the RAFT architecture [34]. We choose this
architecture due to its extensibility to sequential data that we describe in
detail in section 3.3 and visualize some components in figure 2. RAFT uses
three feature embeddings to compute optical flow. First, two feature
embeddings are created with a shared encoder to compute a cost volume, which
is illustrated in blue in figure 2. Second, a context encoder with a different
set of weights creates features, shown in green in figure 2, that are directly
used by a GRU update block [7].

To adopt this idea for event data, we need to build an event representation
that is compatible with CNN layers. Our idea is to split a short sequence of
event data into two subsequences. More specifically, we select an event
sequence ${\cal E}_{i}={\cal E}(t_{i-1},t_{i})$ where $t_{i-1}=t_{i}-\Delta t$
and a second event sequence ${\cal E}_{i+1}={\cal E}(t_{i},t_{i+1})$. For each
event sequence, we create a volumetric voxel grid representation of event data
[41]. It discretizes the time dimension but is able to retain most of the
temporal aspect of events through bilinear interpolation in time. The feature
correlation encoder produces a feature embedding from ${\cal E}_{i}$ and
${\cal E}_{i+1}$ while the context encoder only receives the event
representation of ${\cal E}_{i+1}$. The reason for this is that events between
$t_{i}$ to $t_{i+1}$ contain more accurate contextual information for the
prediction of $F_{i\to i+1}$ than events before time $t_{i}$.

Figure 2: Our proposed recurrent approach for event-based optical flow
estimation. The architecture computes a cost volume $C$ from two subsequent
event representations. Another context encoder $f$ intializes the hidden state
of the GRU update unit and provides additional semantic priors. Our proposed
warping module fuses different timestamps together and facilitates propagation
of information through time.

###  3.3 Differentiable Warm-Starting of Optical Flow

This section elaborates on a method to incorporate flow priors from previous
time steps not only in the evaluation but also in the training of the model.

By default, the RAFT architecture initializes the flow estimate
$F_{i\to{i+1}}$ with 0. Given sequential data, it is also possible to
initialize the optical flow from the last prediction [34]. Different from
prior work, we propose to implement this initialization step in a
differentiable manner and incorporate the warm-starting module in the
training. By doing this, the neural network can incorporate temporal priors
and backpropagate gradients through time. The proposed warm-starting module
does not introduce any additional parameters and can be added to already
trained networks for further finetuning. The overall approach is visualized in
figure 2.

Our goal is to reuse the previous flow prediction $F_{i-1\to i}$ to initialize
the iterative update module that eventually estimates $F_{i\to i+1}$. However,
simply copying the previous flow estimate to initialize the current one would
not take into account that the source pixel location changed according to the
flow itself. Instead, the flow has to be forward warped to the next timestamp.
We achieve this through forward warping of the flow predictions. In forward
warping, multiple source pixels can map to the same target pixel leading to
ambiguities, which is why we perform average splatting in the target pixel
space. We do this by normalizing the warped flow with the weights derived from
the bilinear kernel. In summary, the flow $F_{{i}\to{i+1}}$ is initialized
with the previous flow prediction $F_{{i-1}\to i}$ in the following way:

$\displaystyle g(x_{i})$ $\displaystyle=x_{i}+F_{i\to{i+1}}(x_{i})$ (1)
$\displaystyle F_{i\to{i+1}}(x)$ $\displaystyle=\frac{\sum_{\forall
x_{i-1}}k_{b}(x-g(x_{i-1}))F_{{i-1}\to{i}}(x_{i-1})}{\sum_{\forall
x_{i-1}}k_{b}(x-g(x_{i-1}))}$ (2)

where

$\displaystyle
k_{b}(a)=\text{max}\\{0,1-|a_{x}|\\}\cdot\text{max}\\{0,1-|a_{y}|\\}$ (3)

is the bilinear interpolation kernel. These operations are subdifferentiable
and can directly be implemented in deep learning frameworks such as PyTorch.

###  3.4 Supervision

We supervise all methods on the $L_{1}$ distance between optical flow
prediction and groundtruth. The loss function for our method is defined as

$\displaystyle{\cal
L}=\sum_{i=1}^{N_{i}}\sum_{k=1}^{N_{k}}\gamma^{N_{k}-k}||F_{i\to
i+1}^{\text{gt}}-F_{i\to i+1}^{k}||_{1}$ (4)

where $N_{i}$ is the sequence length, $N_{k}$ the number of iteration steps
per timestep, $F_{i\to i+1}^{\text{gt}}$ the groundtruth for time $i$, and
$F_{i\to i+1}^{k}$ the prediction at the $k$th iteration at time $i$. We set
$\gamma=0.8$.

##  4 DSEC-Flow Dataset

This section introduces our proposed optical flow dataset that overcomes the
major shortcomings of the MVSEC optical flow dataset [40]. Our proposed
dataset extends the publicly available DSEC dataset [11] which contains
disparity ground truth obtained from a rotating Lidar. It features VGA
resolution event cameras which have 3.4 times the number of pixels compared to
the Davis346 used in the MVSEC dataset. Another advantage of the DSEC dataset
is that the disparity maps are cleaned from occlusions, and Lidar points on
moving objects have been removed with a stereo matching filtering technique.
These aspects build a foundation for our extension of this dataset to optical
flow ground truth.

We utilize the intrinsic and extrinsic calibration data for each sequence to
reproject the disparity map into 3D. We then use the odometry ground truth
from the dataset to estimate the forward and backward optical flow at 10 Hz
assuming that the scene is static. However, we noticed that, even though most
lidar points on dynamic objects have been removed via the proposed filtering
technique [11], some points remain on the dynamic objects. To ensure accurate
optical flow ground truth, we manually inspect the whole dataset and extract
all parts of sequences that contain static environments. The resulting dataset
contains 7800 training samples and 2100 test samples in 24 sequences at day
and night. Figure 4 shows examples of optical flow predictions on this
dataset.

An important aspect that indicates the difficulty of an optical flow dataset
is the distribution of pixel displacement magnitudes. Typically, smaller pixel
displacements are easier to predict accurately, in absolute error, than large
displacements [6]. On the other hand, estimating long-range correspondences
has proven to be key for accurate, drift-free optical flow estimation for
high-speed video [17]. Consequently, it is important to devise methods and
datasets that address both small and large displacements.

Figure 3 analyzes displacement distribution of our dataset and compares the
ground truth flow magnitudes of DSEC-Flow with MVSEC. Optical flow ground
truth for MVSEC is originally provided at 20 Hz. However, following initial
evaluation of Zhu et al. [40], recent work also evaluates on temporally
upsampled ground truth at frame rate [19, 41, 37] at approximately 45 Hz.

Figure 3 also shows that our dataset consists of drastically higher
displacement fields than MVSEC. 80 % of optical flow magnitudes on MVSEC are
lower than 4 pixels. For the upsampled version this value decreases to mere
1.7 pixels. In stark contrast, DSEC-Flow reaches 22 pixel flow magnitude for
the 80 % percentile.

For further analysis, we also normalize the flow magnitudes by the image
width. The 80 % percentile for our dataset is 3.4 % while MVSEC reaches 1.3 %
and 0.5 % respectively. This indicates that DSEC-Flow is more challenging even
when normalized with respect to the resolution of the camera.

###  4.1 Evaluation Procedure

Our test set consists of approximately 2100 ground truth flow maps in seven
different sequences. We evaluate optical flow predictions based on the end-
point-error (EPE) and also compute the $N$-pixel errors ($N$PE) that indicate
the percentage of flow errors higher than $N$ pixels in magnitude. We compute
$N$PE for $N$ equals 3, 2 and 1 for coarse to fine analysis of flow outliers,
but consider the EPE as the default metric since it provides the most accurate
measure of flow accuracy.

Figure 3: Cumulative Distribution of flow magnitudes normalized by the number
of pixels in image width. Our dataset consists of significantly higher optical
flow than the MVSEC variants. “ _upsampled_ ” indicates that the dataset
ground truth has been temporally upsampled according to Zhu et. al [40]

##  5 Experiments

We evaluate our method on MVSEC [40] and the proposed dataset DSEC-Flow. The
published baselines are exclusively self-supervised methods that have been
trained and evaluated on the MVSEC dataset. Since we are interested in
understanding the full potential of existing approaches in comparison to our
method, we retrain the EV-FlowNet [40, 41] architecture on both MVSEC and
DSEC-Flow with direct supervision from scratch. Note that the EV-FlowNet
architecture, which was also used in [41], achieved the lowest EPE of all
baselines in the sparse evaluation summarized in table 2.111ECN [37] computes
static optical flow only indirectly via depth and ego-motion. Therefore, the
comparison is not fair. We optimize the learning rate for each method
separately and perform horizonal flipping and random cropping for data
augmention on all datasets. We also set the number of iterations per timestep,
$N_{k}$ in equation 4, to 12 for training and evaluation on both datasets.

###  5.1 MVSEC

|  | Dense |  |  | Sparse |  |
---|---|---|---|---|---|---|---
| Loss Type | EPE | 1PE | 3PE | EPE | 1PE | 3PE
Back to Event Basics [27] | E | - | - | - | 0.92 | - | 5.4
Spike-FlowNet [19] | I | - | - | - | 0.84 | - | -
ECN * [37] | E | 0.35 | - | 0.0 | 0.30 | - | 0.0
EV-FlowNet [40] | I | - | - | - | 0.49 | - | 0.2
Zhu et. al [41] | E | - | - | - | 0.32 | - | 0.0
EV-FlowNet [40] | S | 0.35 | 3.7 | 0.0 | 0.31 | 1.6 | 0.0
Ours | S | 0.27 | 1.7 | 0.0 | 0.24 | 1.7 | 0.0
Table 2: MVSEC results at temporally upsampled optical flow ground truth at
around 45 Hz. Supervision type E refers to event-based loss, I to image-based
(photometric) loss, and S to supervised learning. $*$ indicates that this
method predicts optical flow via depth and ego-motion. Best performance in
bold and second best in italic. We do not highlight the 3PE as it is zero for
almost all methods.

| Dense |  |  | Sparse |  |
---|---|---|---|---|---|---
| EPE | 1PE | 3PE | EPE | 1PE | 3PE
EV-FlowNet [40] | 0.61 | 15.6 | 0.45 | 0.53 | 10.7 | 0.60
Ours | 0.47 | 9.2 | 0.24 | 0.46 | 7.5 | 0.49

Table 3: MVSEC results at original ground truth rate of 20 Hz. Best
performance in bold.

On MVSEC we train our approach with differentiable warm-starting from scratch
with sequence length 2 on the outdoor day 2 sequence. After 10 epochs, we
increase the sequence length to 5. We follow prior work and train, and
evaluate on the temporally upsampled ground truth. The upsampled ground truth
is post-processed to the frame-rate of the camera frames approximately at 45
Hz. Additionally, we also train, and evaluate EV-FlowNet and our method on the
provided ground truth without any modifications to render the task more
challenging.

Note that we also report 3-pixel errors on MVSEC to stay consistent with prior
work in table 2, although 3PE is typically zero.

#### Sparse vs. Dense

Prior work often masks the optical flow predictions to the pixel locations
where events triggered [40, 41, 19, 27]. Since not all methods are open
source, we compare against our approaches and the supervised baseline at
upsampled ground truth. Table 2 shows that the metrics are consistently lower
when the evaluation is performed with masked ground truth. The sparse EPE is
11% lower compared to the dense EPE. Hence, dense optical flow prediction is a
more difficult task on MVSEC. With a dense EPE of 0.27 and sparse EPE of 0.24,
our approach outperforms the sparse and dense baselines and methods proposed
in the literature.

#### Supervised vs. Self-Supervised

Our approach and the retrained EV-FlowNet baseline are the only supervised
networks in the literature on this dataset. Regardeless of sparse or dense
evaluation, supervised methods generally outperform self-supervised methods.
The worse performing supervised baseline, EV-FlowNet, achieves a sparse EPE of
0.31 while Zhu et al. [41] reach 0.32 EPE. We hypothesize that the EPE error
advantage of supervised learning methods would be significantly higher if
MVSEC ground truth optical flow would be of higher quality.

#### Temporally Upsampled vs Original Groundtruth Rate

Table 3 illustrates the performance of supervised EV-FlowNet and our method on
MVSEC without ground truth upsampling. The direct comparison of our method and
the supervised EV-FlowNet on both MVSEC versions reveils that the EPE of our
approach is 23 % lower in both cases. Hence, the original (20 Hz) optical flow
of MVSEC does not yet reveil additional insights with respect to the
performance of estimating larger displacement fields. The experiments of the
next section 5.2 addresses this challenge. Nonetheless, the 1-pixel error
increases substantially for both methods. This is, of course, expected because
the optical flow magnitudes increase due to the larger time windows.

###  5.2 DSEC-Flow

On the DSEC-Flow dataset we perform horizontal flipping and random cropping to
288 pixels in height and 384 pixels in width. We train our non-recurrent
approach with a learning rate of $10^{-4}$ with Adam for 40 epochs until
convergence. Then, we introduce the differentiable warm-starting with a
sequence length of 3 and train for another 6 epochs with random cropping.
Finally, we reduce the learning rate by a factor of 10 and fine-tune on the
full resolution for evaluation. For a fair comparison, we apply the same fine-
tuning schedule to the non-recurrent architecture as well as the baseline
model.

Table 4 compares our approach against the baseline method. Our method achieves
an EPE of 0.79 which is 2.9 times lower than EV-FlowNet with an EPE of 2.32.
To our surprise, EV-FlowNet does not perform comparably in DSEC-Flow. This is
different to the MVSEC experiments where the performance difference between
our method and EV-FlowNet was in a comparable range. A possible explanation is
that EV-FlowNet was designed on MVSEC to work for small pixel displacements.

| EPE | 1PE | 2PE | 3PE
---|---|---|---|---
EV-FlowNet [40] | 2.32 | 55.4 | 29.8 | 18.6
Ours | 0.79 | 12.5 | 4.7 | 2.7
Table 4: Evaluation on the DSEC-Flow dataset. Best performance in bold.

Figure 4: Qualitative examples of optical flow predictions on the DSEC-Flow
test set. Best viewed in PDF format. From left to right: Image for
visualization only, events visualized as image, EV-Flownet, our approach, and
the ground truth. While EV-FlowNet does produce very noisy optical flow, our
proposed approach is able to produce very sharp boundaries consistently.

###  5.3 Ablation Study

This section introduces the ablation studies that highlight the significance
of different modules. We perform these ablation studies both on the DSEC-Flow
dataset and on the MVSEC dataset. Table 5 summarizes the results of the
ablation study that focuses on warm-starting strategies in combination with
different evaluation and training combinations.

| DSEC-Flow | MVSEC
---|---|---
| EPE | 1PE | 2PE | 3PE | EPE | 1PE
Ours w/o WS | 0.82 | 13.5 | 5.2 | 3.0 | 0.53 | 12.2
Ours eval WS | 0.88 | 14.1 | 5.6 | 3.3 | 0.51 | 10.9
Ours w/o WS: 100 iters | 0.87 | 14.9 | 5.9 | 3.4 | 0.60 | 16.3
Ours train & eval WS: 100 iters | 0.88 | 15.1 | 6.0 | 3.5 | 0.49 | 10.0
Ours train & eval WS | 0.79 | 12.7 | 4.7 | 2.7 | 0.47 | 9.2
Table 5: Ablation study on DSEC-Flow and MVSEC. Best performance in bold. _WS_
refers to warm-starting optical flow from the previous estimated optical flow.
_eval WS_ means that warm-starting is only used at inference time without
training with warm-starting beforehand.

The baseline method (Ours w/o WS) is the proposed non-recurrent architecture
applied to event data and does not incorporate any warm-starting in training
or evaluation.

#### DSEC-Flow Dataset

On the DSEC-Flow dataset, the baseline approach performs better compared to
evaluating the same method with warm-starting. That is without incorporating
warm-starting in the training. Surprisingly, it also outperforms both
evaluation protocols that involve 100 iterations per timestep, regardless
whether we train and evaluate with warm-starting or not. Our proposed
approach, which incorporates warm-starting in the training and evaluation,
outperforms the ablated methods with 4.0 % lower EPE of 0.79.

#### MVSEC Dataset

The results on the MVSEC ablation similarly suggest that training and
evaluating with warm-starting is superior to the aforementioned variations.
However, differently from DSEC-Flow, purely evaluating with warm-starting
(without training) leads to lower errors than not using warm-starting at all.
Still, the model that also incorporates training with the warm-starting module
achieves 4.2 % lower EPE than the closest ablated method (training +
evaluation + 100 iterations).

Overall, the effect of incorporating warm-starting into the fine-tuning
procedure improves the performance consistently by a few percent on both
datasets. Surprisingly, such gains are not consistent in our experiments if we
only add warm-starting at inference time as reported in [34].

##  6 Conclusion

This work proposes an alternative approach to event-based optical flow
estimation by adopting design choices from frame-based optical flow estimation
and modifying the method to integrate with event-based input. We showed that
the proposed method drastically outperforms baseline methods on two datasets,
which suggests that the commonly used U-Net architecture in the event-based
literature is a suboptimal choice for the task at hand. A key contribution of
this work is a novel real-world optical flow dataset that addresses the major
shortcomings of existing event-based datasets. We showed that training on this
dataset leads to high-quality optical flow predictions that are unattainable
with previous event-based vision datasets.

##  7 Acknowledgments

This work was supported by the National Centre of Competence in Research
(NCCR) Robotics through the Swiss National Science Foundation (SNSF) and by
the European Research Council (ERC) under Grant Agreement 864042 (AGILEFLIGHT)

## References

  * [1] Mohammed Almatrafi, Raymond Baldwin, Kiyoharu Aizawa, and Keigo Hirakawa.  Distance surface for event-based optical flow.  IEEE transactions on pattern analysis and machine intelligence, 2020\.
  * [2] Patrick Bardow, Andrew J. Davison, and Stefan Leutenegger.  Simultaneous optical flow and intensity estimation from an event camera.  In IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), pages 884–892, 2016.
  * [3] Ryad Benosman, Charles Clercq, Xavier Lagorce, Sio-Hoi Ieng, and Chiara Bartolozzi.  Event-based visual flow.  IEEE Trans. Neural Netw. Learn. Syst., 25(2):407–417, 2014.
  * [4] Ryad Benosman, Sio-Hoi Ieng, Charles Clercq, Chiara Bartolozzi, and Mandyam Srinivasan.  Asynchronous frameless event-based optical flow.  Neural Networks, 27:32–37, 2012.
  * [5] Tobias Brosch, Stephan Tschechne, and Heiko Neumann.  On event-based optical flow detection.  Frontiers in neuroscience, 9:137, 2015.
  * [6] D. J. Butler, J. Wulff, G. B. Stanley, and M. J. Black.  A naturalistic open source movie for optical flow evaluation.  In A. Fitzgibbon et al. (Eds.), editor, European Conf. on Computer Vision (ECCV), Part IV, LNCS 7577, pages 611–625. Springer-Verlag, Oct. 2012.
  * [7] Kyunghyun Cho, B van Merrienboer, Caglar Gulcehre, F Bougares, H Schwenk, and Yoshua Bengio.  Learning phrase representations using rnn encoder-decoder for statistical machine translation.  In Conference on Empirical Methods in Natural Language Processing (EMNLP 2014), 2014.
  * [8] Alexey Dosovitskiy, Philipp Fischer, Eddy Ilg, Philip Häusser, Caner Hazırbaş, Vladimir Golkov, Patrick van der Smagt, Daniel Cremers, and Thomas Brox.  FlowNet: Learning optical flow with convolutional networks.  In Int. Conf. Comput. Vis. (ICCV), pages 2758–2766, 2015.
  * [9] Guillermo Gallego, Tobi Delbruck, Garrick Orchard, Chiara Bartolozzi, Brian Taba, Andrea Censi, Stefan Leutenegger, Andrew Davison, Jörg Conradt, Kostas Daniilidis, and Davide Scaramuzza.  Event-based vision: A survey.  IEEE Trans. Pattern Anal. Mach. Intell., 2020.
  * [10] Daniel Gehrig, Henri Rebecq, Guillermo Gallego, and Davide Scaramuzza.  EKLT: Asynchronous photometric feature tracking using events and frames.  Int. J. Comput. Vis., 2019.
  * [11] Mathias Gehrig, Willem Aarents, Daniel Gehrig, and Davide Scaramuzza.  Dsec: A stereo event camera dataset for driving scenarios.  IEEE Robot. Autom. Lett., 2021.
  * [12] Andreas Geiger, Philip Lenz, Christoph Stiller, and Raquel Urtasun.  Vision meets robotics: The KITTI dataset.  Int. J. Robot. Research, 32(11):1231–1237, 2013.
  * [13] Berthold K.P. Horn and Brian G. Schunck.  Determining optical flow.  J. Artificial Intell., 17(1):185 – 203, 1981.
  * [14] Tak-Wai Hui, Xiaoou Tang, and Chen Change Loy.  Liteflownet: A lightweight convolutional neural network for optical flow estimation.  In IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), 2018.
  * [15] Junhwa Hur and Stefan Roth.  Iterative residual refinement for joint optical flow and occlusion estimation.  In IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), 2019.
  * [16] Joel Janai, Fatma Güney, Aseem Behl, Andreas Geiger, et al.  Computer vision for autonomous vehicles: Problems, datasets and state of the art.  Foundations and Trends® in Computer Graphics and Vision, 12(1–3):1–308, 2020.
  * [17] Joel Janai, Fatma Güney, Jonas Wulff, Michael Black, and Andreas Geiger.  Slow flow: Exploiting high-speed cameras for accurate and diverse optical flow reference data.  In CVPR, 2017.
  * [18] Nima Khademi Kalantari and Ravi Ramamoorthi.  Deep high dynamic range imaging of dynamic scenes.  ACM Trans. Graph., 36(4):144–1, 2017.
  * [19] Chankyu Lee, Adarsh Kumar Kosta, Alex Zihao Zhu, Kenneth Chaney, Kostas Daniilidis, and Kaushik Roy.  Spike-flownet: event-based optical flow estimation with energy-efficient hybrid neural networks.  In Eur. Conf. Comput. Vis. (ECCV), 2020.
  * [20] Patrick Lichtsteiner, Christoph Posch, and Tobi Delbruck.  A 128$\times$128 120 dB 15 $\mu$s latency asynchronous temporal contrast vision sensor.  IEEE J. Solid-State Circuits, 43(2):566–576, 2008.
  * [21] Pengpeng Liu, Michael Lyu, Irwin King, and Jia Xu.  Selflow: Self-supervised learning of optical flow.  In IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), 2019.
  * [22] Bruce D. Lucas and Takeo Kanade.  An iterative image registration technique with an application to stereo vision.  In Int. Joint Conf. Artificial Intell. (IJCAI), pages 674–679, 1981\.
  * [23] Moritz Menze and Andreas Geiger.  Object scene flow for autonomous vehicles.  In IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), 2015.
  * [24] Elias Mueggler, Christian Forster, Nathan Baumli, Guillermo Gallego, and Davide Scaramuzza.  Lifetime estimation of events from dynamic vision sensors.  In IEEE Int. Conf. Robot. Autom. (ICRA), pages 4874–4881, 2015\.
  * [25] Simon Niklaus and Feng Liu.  Softmax splatting for video frame interpolation.  In IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), June 2020.
  * [26] Liyuan Pan, Miaomiao Liu, and Richard Hartley.  Single image optical flow estimation with an event camera.  IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), 2020.
  * [27] Federico Paredes-Vallés and Guido CHE de Croon.  Back to event basics: Self-supervised learning of image reconstruction for event cameras via photometric constancy.  arXiv preprint arXiv:2009.08283, 2020.
  * [28] Tong Qin, Peiliang Li, and Shaojie Shen.  Vins-mono: A robust and versatile monocular visual-inertial state estimator.  IEEE Transactions on Robotics, 34(4):1004–1020, 2018.
  * [29] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [30] Bodo Rueckauer and Tobi Delbruck.  Evaluation of event-based algorithms for optical flow with ground-truth from inertial measurement sensor.  Frontiers in neuroscience, 2016.
  * [31] Karen Simonyan and Andrew Zisserman.  Two-stream convolutional networks for action recognition in videos.  In Conf. Neural Inf. Process. Syst. (NeurIPS), 2014.
  * [32] Deqing Sun, Xiaodong Yang, Ming-Yu Liu, and Jan Kautz.  PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume.  In IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), 2018.
  * [33] Deqing Sun, Xiaodong Yang, Ming-Yu Liu, and Jan Kautz.  Models matter, so does training: An empirical study of cnns for optical flow estimation.  IEEE transactions on pattern analysis and machine intelligence, 2019\.
  * [34] Zachary Teed and Jia Deng.  Raft: Recurrent all-pairs field transforms for optical flow.  In Eur. Conf. Comput. Vis. (ECCV), 2020.
  * [35] Jia Xu, Rene Ranftl, and Vladlen Koltun.  Accurate optical flow via direct cost volume processing.  In IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), July 2017.
  * [36] Gengshan Yang and Deva Ramanan.  Volumetric correspondence networks for optical flow.  NeurIPS, 5:12, 2019.
  * [37] C. Ye, A. Mitrokhin, C. Fermüller, J. A. Yorke, and Y. Aloimonos.  Unsupervised learning of dense optical flow, depth and egomotion with event-based sensors.  In 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2020.
  * [38] Chengxi Ye, Anton Mitrokhin, Chethan Parameshwara, Cornelia Fermüller, James A. Yorke, and Yiannis Aloimonos.  Unsupervised learning of dense optical flow and depth from sparse event data.  arXiv e-prints, 2019.
  * [39] Shengyu Zhao, Yilun Sheng, Yue Dong, Eric I Chang, Yan Xu, et al.  Maskflownet: Asymmetric feature matching with learnable occlusion mask.  In IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), 2020.
  * [40] Alex Zihao Zhu, Liangzhe Yuan, Kenneth Chaney, and Kostas Daniilidis.  EV-FlowNet: Self-supervised optical flow estimation for event-based cameras.  In Robotics: Science and Systems (RSS), 2018.
  * [41] Alex Zihao Zhu, Liangzhe Yuan, Kenneth Chaney, and Kostas Daniilidis.  Unsupervised event-based learning of optical flow, depth, and egomotion.  In IEEE Conf. Comput. Vis. Pattern Recog. (CVPR), 2019.


aainstitutetext: CAS Key Laboratory of Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, P. R. China
bbinstitutetext: School of Physical Sciences, University of Chinese Academy of
Sciences, Beijing 100049, P. R. Chinaccinstitutetext: Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing 100049, P. R. China

# Enhanced Higgs pair production from higgsino decay at the HL-LHC

Jianpeng Dai c,b Tao Liu a,b Daohan Wang a,b Jin Min Yang
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The scenario of multi-sector SUSY breaking predicts pseudo-goldstinos which
are not absorbed by the gravitino and their mass can be as low as ${{\cal
O}(0.1)}$ GeV. Since the interactions of pseudo-goldstinos are not so weak as
gravitino, a produced higgsino can decay to a pseudo-goldstino plus a Higgs
boson insider the detector at the LHC, and thus the higgsino pair production
can lead to the signal of Higgs pair plus missing energy. For the scenario of
natural SUSY which requires rather light higgsinos, such events may sizably
outnumber the Higgs pair events predicted by the SM and be accessible at the
HL-LHC (14 TeV with a luminosity of 3 $\rm{ab}^{-1}$). In this work we examine
the observability of such Higgs pair plus missing energy from the decay of
light higgsinos produced at the HL-LHC. Considering three channels of the
Higgs-pair decay ($bbWW^{*}$, $bb\gamma\gamma$, $bbbb$), our detailed Monte
Carlo simulations for the signal and backgrounds show that the best channel is
$bbbb+{\>/\\!\\!\\!\\!E_{T}}$, whose statistical significance can reach
$2\sigma$ level for a light higgsino allowed by current experiments. This is
over the SM Higgs pair result which is about $1.8\sigma$.

##  1 Introduction

The precision test of the Higgs properties and the exploration of new physics
beyond the standard model (SM) is the main task of the Large Ladron Collider
(LHC) and other future colliders after the Higgs boson was discovered in 2012.
Now the Higgs mass has already been measured to an impressive precision, and
in order to confirm the mechanism of electroweak symmetry breaking described
in the SM and probe possible new physics beyond the SM, the Higgs self-
coupling $\lambda_{hhh}$ has to be reconstructed by experimental measurements
which requires the production of at least two Higgs bosons at the LHC. In the
SM, the dominant channel for Higgs pair production at the LHC is the gluon
fusion channel, in which these two Higgs bosons could be directly radiated
from heavy quark loops or from a virtual Higgs boson through the self-coupling
$\lambda_{hhh}$. Comparing with the single Higgs production, the Higgs pair
processes are additionally suppressed by the destructive interference between
the above two contributions. The cross sections for Higgs pair production is
about three orders of magnitude smaller than the single Higgs production, and
thus it cannot reach the observable level when considering the huge QCD
background at the LHC. On the experimental side, ATLAS and CMS have published
measurements of Higgs pair production in the decay channels $b{\bar{b}}WW^{*}$
Aaboud:2018zhh ; CMS:2017ums , $b{\bar{b}}b{\bar{b}}$ Aaboud:2018knk ;
CMS:2018smw , $b{\bar{b}}\gamma\gamma$ Sirunyan:2018iwt ; Aaboud:2018ftw ,
$b{\bar{b}}\tau^{+}\tau^{-}$ Aaboud:2018sfw ; Sirunyan:2017djm  and
$WW^{*}WW^{*}$ Aaboud:2018ksn , and the combination of different measurements
gives a constraint on the triliear self-coupling,
$-5.0<\lambda_{hhh}/\lambda_{hhh}^{SM}<12.0$ at $95\%$ C.L. 2020cos , with the
assumption that all other couplings are SM-like. The high-luminosity run of
the LHC will further narrow the window of $\lambda_{hhh}$.

The Higgs pair production rate could be altered in new physics models (see,
for examples Refs. Contino:2010mh ; Grober:2010yv ; Dolan:2012ac ; Cao:2013si
; Han:2013sga ; Hespel:2014sla ; Cao:2014kya ; Dawson:2015oha ; Lu:2015qqa ;
Kang:2015nga ; 2016yan ; 2017yan ; 2017caoyan ; Nakamura:2017irk ;
Chang:2017niy ; Ren:2017jbg ; Huang:2017nnw ; Adhikary:2017jtu ;
Agrawal:2019bpm ; Cheung:2020xij ; 2020yan ; Abouabid:2021yvw  and references
therein) which may predict different Higgs self-coupling from the SM, e.g., in
the next-to-minimal supersymmetric model (NMSSM) or little Higgs theory, the
Higgs self-coupling may be significantly different from the SM value
Wu:2015nba ; Wang:2007zx ; Han:2009zp . Of course, the Higgs pair events may
also come from the decays of the new particles predicted in new physics
models. In this work we consider supersymmetric models with multi-sector SUSY
breaking, in which each breaking sector provides a massless goldstino
$\eta_{i}$ at tree level with SUSY breaking scale $F_{i}$. When taking into
account radiative corrections, the true goldstino which will be absorbed by
gravitino remains massless, while the pseudo-goldstino acquires non-vanishing
contribution and obtains a mass Argurio:2011hs ; Dai:2021eah . Unconstrained
by the supercurrent like gravitino, the interaction between pseudo-goldstino
and visible fields could be strong enough and thus may lead to unconventional
phenomenology Cheung:2010mc ; Cheung:2010qf ; Craig:2010yf ; Thaler:2011me ;
Cheung:2011jq ; Argurio:2011gu ; Argurio:2011hs ; Dudas:2011kt ; Liu:2013sx ;
Hikasa:2014yra ; Liu:2014lda ; Ferretti:2013wya ; Dai:2021eah ;
Franzosi:2021zwp . In our previous studies Hikasa:2014yra ; Liu:2014lda ;
Liu:2013sx  we found that the lightest neutralino could be bino-like and
mainly decay to Higgs or longitudinal $Z$ boson plus pseudo-goldstino while
the lightest chargino may decay to $W$ boson plus pseudo-goldstino. In those
studies only the decay channels to $Z$ or $W$ boson are considered since they
are easier to detect. In case that the lightest neutralino dominantly decays
to Higgs plus pseudo-goldstino, the production of neutralino pair will lead to
the signal of Higgs pair plus missing energy. Note that as fermions the
neutralinos may be more copiously produced than the Higgs bosons. The Higgs
pair events from these neutralino decays may sizably outnumber the Higgs pair
events predicted by the SM. In addition, the accompanied large missing energy
may help to reduce the QCD backgrounds. In this work we perform a
comprehensive study for this Higgs channel at the high luminosity LHC (HL-
LHC). We will work in the natural SUSY scenario as an example which predicts
light higgsinos, and consider three typical Higgs-pair decay channels to study
the observability, i.e., $b{\bar{b}}WW^{*}$ with $W$ boson decaying
leptonically, $b{\bar{b}}b{\bar{b}}$ and $b{\bar{b}}\gamma\gamma$.

This work is organized as follows. In Sect. 2 we make a brief review on
pseudo-goldstino and mention some related theoretical basics. Then in Sect. 3
we perform a Monte Carlo simulation for the signal and backgrounds for three
decay channels for the Higgs pair. The conclusion is made in Sect. 4.

##  2 A description of theoretical framework

In the scenario of multi-sector SUSY breaking, each hidden sector with
spontaneous SUSY breaking at scale $F_{i}$ could be parameterized in a non-
linear way
$X_{i}=\eta_{i}^{2}/(2F_{i})+\sqrt{2}\theta\eta_{i}+\theta^{2}F_{i}$, where
$\eta_{i}$ denotes the so-called goldstino. Then the soft terms for the
visible superfields can be obtained through the non-trivial Kähler potential
$K$ and gauge kinetic function $f$

$\displaystyle K$ $\displaystyle=$
$\displaystyle\Phi^{\dagger}\Phi\sum_{i}\frac{m_{\phi,i}^{2}}{F_{i}^{2}}X_{i}^{\dagger}X_{i},$
(1) $\displaystyle f_{ab}$ $\displaystyle=$
$\displaystyle\frac{1}{g_{a}^{2}}\delta_{ab}\left(1+\sum_{i}\frac{2m_{a,i}}{F_{i}}X_{i}\right).$
(2)

Here $m_{\phi,a}$ are the soft masses for the chiral superfields and gauginos,
respectively. Other soft trilinear $A$ terms and bilinear $B_{\mu}$ which will
not be used in the following discussion could also be easily constructed. Take
the scenario of two hidden sectors as an example, with the definition
$F=\sqrt{F_{1}^{2}+F_{2}^{2}}$ and $\tan\theta=F_{2}/F_{1}$,
$G=\eta_{1}\cos\theta+\eta_{2}\sin\theta$ would be absorbed by gravitino
through the super-Higgs mechanism and
$G^{\prime}=-\eta_{1}\sin\theta+\eta_{2}\cos\theta$ is left as the physical
pseudo-goldstino. Then one can obtain the interaction between $G,G^{\prime}$
and visible particles up to order $1/F_{i}$

$\displaystyle\mathcal{L}_{G}$ $\displaystyle=$
$\displaystyle\frac{m_{\phi}^{2}}{F}G\psi\phi^{*}-\frac{im_{a}}{\sqrt{2}F}G\sigma^{\mu\nu}\lambda^{a}F^{a}_{\mu\nu}+\frac{m_{a}}{F}G\lambda^{a}D^{a},$
(3) $\displaystyle\mathcal{L}_{G^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{\widetilde{m}_{\phi}^{2}}{F}G^{\prime}\psi\phi^{*}-\frac{i\widetilde{m}_{a}}{\sqrt{2}F}G^{\prime}\sigma^{\mu\nu}\lambda^{a}F^{a}_{\mu\nu}+\frac{\widetilde{m}_{a}}{F}G^{\prime}\lambda^{a}D^{a}.$
(4)

Through the definition of $m_{a,\phi}$ and $\widetilde{m}_{a,\phi}$

$\displaystyle
m_{a}=m_{a,1}+m_{a,2},\quad\widetilde{m}_{a}=-m_{a,1}\tan\theta+m_{a,2}\cot\theta,$
$\displaystyle
m^{2}_{\phi}=m^{2}_{\phi,1}+m^{2}_{\phi,2},\quad\widetilde{m}^{2}_{\phi}=-m^{2}_{\phi,1}\tan\theta+m^{2}_{\phi,2}\cot\theta,$
(5)

it is easy to see that pseudo-goldstino $G^{\prime}$ could couple to ordinary
fields in a total different way compared with gravitino. Under the condition
of approximately vanishing $\widetilde{m}_{a}$, the lightest neutralino could
only decay to Higgs or longitudinal $Z$ boson plus $G^{\prime}$. The details
about the derivation as well as the corresponding dynamical realization of
this condition from the viewpoint of model-buildings could be found in our
previous paper Hikasa:2014yra .

Obviously the mass of $G^{\prime}$ is a crucial parameter for phenomenological
analysis. At tree level $G^{\prime}$ acquires contribution which is just twice
the gravitino mass $m_{3/2}$ via supergravity Cheung:2010mc . The leading
order correction without gravitational effect arise at three-loop level. When
the two messenger scales are equal, the loop correction is at the order of GeV
scale Argurio:2011hs . And through an explicit analytical calculation, we
found that this correction could be as low as ${\cal O}(0.1)$ GeV if one
messenger mass is 100 times larger than the other Dai:2021eah . In our
following Motel Carlo simulations the pseudo-goldstino mass is fixed to a
reasonable value of $0.5$ GeV.

Next, let us turn to the visible sector and have a close look at neutralino
pair production. The property of the neutralinos depends heavily on the choice
of certain parameters in the minimal supersymmetric model (MSSM). Since our
aim is to study the Higgs phenomenology arising from pseudo-glodstino, we
choose the scenario of natural SUSY Brust:2011tb ; Papucci:2011wy ;
Hall:2011aa ; Feng:2012jfa ; Baer:2012up ; Tata:2020afe  which predicts light
higgsinos. In this scenario the higgsinos are rather light, i.e., $\mu$ is
assumed to be 100-300 GeV, while the gauginos are much heavier. So the
neutralinos $\tilde{\chi}^{0}_{1}$ and $\tilde{\chi}^{0}_{2}$ are highly
higgsino-like and nearly degenerate in masses. Then at the LHC the productions
of any pair of them just give missing energy in the MSSM. In order to detect
their productions, a hard jet radiated from initial partons is usually
required and the signal of monojet plus missing energy is searched, which is
found to be rather challenging at the LHC Han:2013usa  (a global likelihood
analysis of the electroweakino sector shows that no range of neutralino or
chargino masses can be robustly excluded by current LHC searches
GAMBIT:2018gjo ). Now with multi-sector SUSY breaking, the higgsino-like
$\tilde{\chi}^{0}_{1}$ and $\tilde{\chi}^{0}_{2}$ could decay to Higgs boson
plus pseudo-goldstino, and their pair productions can lead to the signal of
Higgs pair plus missing energy.

Figure 1:  The Feynman diagram for the Higgs pair event from the decay of
higgsinos produced at the LHC

For the pair productions of $\tilde{\chi}^{0}_{1}$ and $\tilde{\chi}^{0}_{2}$
at the LHC, the cross-section of $\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{2}$
production is much larger than $\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}$ and
$\tilde{\chi}^{0}_{2}\tilde{\chi}^{0}_{2}$. Thus

$\displaystyle
pp\rightarrow\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{2}\rightarrow
hhG^{\prime}G^{\prime}$ (6)

is considered in our study. The corresponding Feynman diagram of the signal is
shown in Fig. 1. To get the mass spectrum and the corresponding mixing
matrices for the neutralinos, SOFTSUSY 2002  is used in our calculation. Here
the mass of the SM-liked Higgs boson is fixed at $125$ GeV, $\mu$ is assumed
to have a positive sign and the value of $\tan\beta$ is chosen to be 10. In
addition, we fix soft gaugino masses as $M_{1}=1$ TeV and $M_{2}=1$ TeV.

Figure 2:  The Higgs pair cross section from
$\tilde{\chi}_{1}^{0}\tilde{\chi}^{0}_{2}$ production at the 14 TeV LHC at
next-to-leading order, compared with the SM Higgs pair cross-section at NLO
Borowka:2016ehy ; Baglio:2018lrj .

Before discussing the role of missing energy from pseudo-goldstino in
distinguishing signal from backgrounds, in Fig. 2 we show the total cross
sections for Higgs pair produced through higgsino decay in natural SUSY,
compared with the SM Higgs pair cross-section at NLO Borowka:2016ehy ;
Baglio:2018lrj . For the couplings between higgsinos and pseudo-goldstino, we
use an effective way to study the phenomenology as in the literature
Hikasa:2014yra  and assume that the branching ratio of higgsino decay to Higgs
to be 1/2 throughout this paper. So we see that for light higgsinos the Higgs
pair produced from higgsino decay may have a larger cross section than the SM
Higgs pair.

##  3 Observability at the HL-LHC

In this section we study in detail the signal of Higgs pair plus missing
energy from higgsino decay in natural SUSY. Signal and the background
processes are modeled using simulated Monte Carlo event samples by MadGraph5_
<EMAIL_ADDRESS>Alwall:2014hca  with the default NNPDF2.3QED parton distribution
functions Ball:2013hta  at the $\sqrt{s}=14$ TeV LHC. The cross section of
higgsino production is normalized to NLO with the help of Prospino2
Beenakker:1999xh . The pseudo-goldstino interaction is implemented in
FeynRules 2014FR  and the UFO model file 2012UF  is passed to MadGraph5. We
use PYTHIA8.205 Sjostrand:2014zea  program to describe the parton-shower and
hadronization.

The fast detector simulations are performed by Delphes deFavereau:2013fsa
with the ATLAS detector. Using FastJet Cacciari:2011ma  for jet-reconstruction
with the anti-$k_{T}$ algorithm Cacciari:2008gp , fixing a cone size of $R$=
0.4 for a jet (the situation of $bbbb$ channel will be further discussed in
the following subsection), we include detector effects relevant for the HL-
LHC, where jets and leptons are smeared according to their energies. For the
analysis, we consider jets with $p_{Tj}>20$ GeV and $|\eta_{j}|<2.5$, and use
the flat $b$-tagging efficiency $\epsilon_{b\rightarrow b}$= 0.75 and the flat
mis-tagging rates for non-$b$ jets $\epsilon_{c\rightarrow b}$= 0.1 and
$\epsilon_{j\rightarrow b}$= 0.01 ATL-PHYS-PUB-2019-005 . As for leptons, we
require $P^{\ell}_{T}/\sum P_{T}>0.7$ with $P_{T}^{\ell}>20$ GeV and
$|\eta_{\ell}|<2.5$ within $\Delta R$=0.3.

###  3.1 The signal of $hhG^{\prime}G^{\prime}\rightarrow
b\bar{b}WW^{*}G^{\prime}G^{\prime}\rightarrow
b\bar{b}\ell^{+}\ell^{-}+{\>/\\!\\!\\!\\!E_{T}}$

We first study the channel $hhG^{\prime}G^{\prime}\rightarrow
b\bar{b}WW^{*}G^{\prime}G^{\prime}\rightarrow
bb\ell^{+}\ell^{-}+{\>/\\!\\!\\!\\!E_{T}}$, which have exactly the same
signature as the SM higgs pair production channel $hh\rightarrow
b\bar{b}WW^{*}\rightarrow bb\ell^{+}\ell^{-}+{\>/\\!\\!\\!\\!E_{T}}$. The
branching ratio of $h\rightarrow WW^{*}$ is the second largest, next to
$h\rightarrow b\bar{b}$. Considering all relevant branching fractions, our
signal cross section is
$\frac{1}{2}\cdot\sigma_{\tilde{\chi}_{1}^{0}\tilde{\chi}_{2}^{0}}\cdot
BR(h\rightarrow WW^{*}\rightarrow\ell^{+}\ell^{-}\nu\bar{\nu})\cdot
BR(h\rightarrow b\bar{b})$, where $\ell$ denotes an electron or a muon,
including leptons from tau decay. Because our signal is similar to SM Higgs
pair channels, so we can refer to their backgrounds directly 2019PDH . The
major background is $t\bar{t}$ production, whose NNLO QCD cross-section is
953.6 pb Czakon:2013goa . The second large background is $t\bar{t}h$, whose
NLO QCD cross-section is 611.3 fb Dittmaier:2011ti . For the
$t\bar{t}V(V=W^{\pm},Z)$ backgrounds, we apply an NLO k-factor of 1.54 and
obtain a cross-section of 1.71 pb deFlorian:2016spz . We also apply a k-factor
of 1.0 for the Drell-Yan type background $\tau\tau bb$. Finally, we generate
$tW^{\pm}j$ events (in the five flavor scheme), whose overlap with $t\bar{t}$
should be subtracted. To reconstruct events, the off-shell effects for the top
quark and $W$ boson need to be considered properly.

We employ the following cuts at parton level: $p_{Tj}>20$ GeV, $p_{Tb}>20$
GeV, $p_{T\gamma}>10$ GeV, $p_{T\ell}>10$ GeV, $|\eta_{j}|$ < 5, $|\eta_{b}|$
< 5, $|\eta_{\gamma}|$ < 2.5, $|\eta_{\ell}|$ < 2.5, $\Delta R_{bb}<$ 1.8,
$\Delta R_{\ell\ell}<$ 1.3, 70 GeV $<m_{jj},m_{bb}<$ 160 GeV and
$m_{\ell\ell}<$ 75 GeV. For $tW^{\pm}j$ backgrounds, we impose 5 GeV
$<m_{\ell\ell}<$ 75 GeV additionally. Here the angular distance $\Delta
R_{ij}$ is defined by

$\displaystyle\Delta
R_{ij}=\sqrt{(\Delta\phi_{ij})^{2}+(\Delta\eta_{ij})^{2}},$ (7)

where $\Delta\phi_{ij}=\phi_{i}-\phi_{j}$ and
$\Delta\eta_{ij}=\eta_{i}-\eta_{j}$ are respectively the differences of the
azimuthal angles and rapidities between particles $i$ and $j$.

Figure 3:  Distribution of the $\mu$ =150 GeV for signal and different types
of backgrounds after basic cuts. The $y$-axis represents the normalized number
of events for each process.

We set a sequence of event selections. From Fig. 3 we can find that the
distributions of $\Delta R_{\ell\ell}$, $\Delta R_{bb}$, $m_{bb}$,
$m_{\ell\ell}$ and ${\>/\\!\\!\\!\\!E_{T}}$ are different between the signal
and the backgrounds, and thus we use them to make further cuts. Since the
$P_{T}$ distributions of backgrounds and signal have no obvious difference, we
simply set $P_{Tl}>20$ GeV and $P_{Tb}>30$ GeV. The cuts $\Delta R_{ll}<1.0$
and $\Delta R_{bb}<1.0$ are imposed to highly suppress the backgrounds. We
make use of invariant masses of leptons and $b-$jets to suppress backgrounds,
because the pair of $b-$jets is reconstructed near the Higgs boson mass. From
a further analysis of the signal and background features, we find that the
missing energy of $t\bar{t}$ is quite approaching to the signal, and thus we
gradually increase the cuts of ${\>/\\!\\!\\!\\!E_{T}}$ to find the best
value. Finally, we find that ${\>/\\!\\!\\!\\!E_{T}}$ > 120 GeV can most
efficiently cut the $t\bar{t}$ and $t\bar{t}h$ backgrounds. The detailed
results are shown in Table 1.

Table 1:  Signal and background cross sections in units of fb after baseline
cuts (first row) and at different stages of analysis, using a combination of
kinematic variables and requiring the events number $N>10$. The significance
$\sigma$ is calculated using the Poisson formula for a luminosity of 3
$\rm{ab}^{-1}$ at the 14 TeV LHC.

Cuts | signal | $t\bar{t}$ | $t\bar{t}h$ | $t\bar{t}V$ | $tW+j$ | $b\bar{b}\tau^{+}\tau^{-}$ | $\sigma$ | $S/B$
---|---|---|---|---|---|---|---|---
$P_{T\ell}$ > 20 GeV, $P_{Tb}$ > 30 GeV | 0.1008 | 3593.325 | 5.0506 | 10.684 | 144.19 | 1.169 | 0.01 | 2.7$\times 10^{-6}$
$\Delta R_{\ell\ell}<1.0,~{}~{}\Delta R_{bb}<1.0$ | 0.015 | 30.59 | 0.065 | 0.175 | 0.9424 | 0.0254 | 0.146 | 4.7$\times 10^{-4}$
$m_{\ell\ell}$ < 65 GeV, 95 GeV $<m_{bb}<140$ GeV | 0.00784 | 0.4746 | 0.0165 | 0.01287 | 0.0228 | 0.0039 | 0.58 | 0.0148
${\>/\\!\\!\\!\\!E_{T}}$ > 50 GeV | 0.00682 | 0.4 | 0.0144 | 0.0116 | 0.023 | 0.0039 | 0.553 | 0.015
${\>/\\!\\!\\!\\!E_{T}}$ > 80 GeV | 0.0058 | 0.268 | 0.01 | 0.0099 | 0.0228 | 0.0039 | 0.564 | 0.018
${\>/\\!\\!\\!\\!E_{T}}$ > 120 GeV | 0.00433 | 0.0603 | 0.006 | 0.0064 | 0.0114 | 0.0039 | 0.793 | 0.049

With the above analysis, we can summarize the sequence of cuts as follows:

  * •

The two leading jets must be $b-$tagged, each with $P_{T}$ > 30 GeV;

  * •

Exactly two isolated leptons of opposite sign, each with $P_{T\ell}$ > 20 GeV;

  * •

Proximity cut of $\Delta R{{}_{\ell\ell}}$ < 1.0 for the two leptons;

  * •

Proximaty cut of $\Delta R{{}_{bb}}$ < 1.0 for two $b-$tagged jets;

  * •

$m_{\ell\ell}$ < 65 GeV for two leptons;

  * •

95 GeV < $m_{bb}$ < 140 GeV for the two $b-$tagged jets;

  * •

${\>/\\!\\!\\!\\!E_{T}}=|{\;/\\!\\!\\!\\!\vec{P}_{T}}|>120$ GeV for the
reconstructed missing transverse momentum.

From our simulation result, we find the $bb\tau\tau$ background can be greatly
suppressed, while the $t\bar{t}$ background remains to be the dominant one.
Finally the significance can only reach 0.79 $\sigma$. Although we can
increase the value of $\mu$ to relatively enhance the signal events with
${\>/\\!\\!\\!\\!E_{T}}>120$ GeV, this will also suppress the total signal
cross section and thus cannot enhance the significance. So we conclude that
this channel cannot be observed at the HL-LHC.

###  3.2 The signal of $hhG^{\prime}G^{\prime}\rightarrow
bb\bar{b}\bar{b}+{\>/\\!\\!\\!\\!E_{T}}$

Now we turn to the decay channel $hhG^{\prime}G^{\prime}\rightarrow
bb\bar{b}\bar{b}+{\>/\\!\\!\\!\\!E_{T}}$. The decay $h\rightarrow b\bar{b}$
has the largest branching ratio of $58\%$. In the SM this channel from the
Higgs pair production suffers from overwhelming QCD multi-jet backgrounds and
its significance can only reach $1.8\sigma$ at the HL-LHC 2015NHP . In the
natural SUSY scenario under our consideration, ${\>/\\!\\!\\!\\!E_{T}}$ of the
signal has the same order of magnitude as $\mu$, so the $bbbb$ process could
be one of the main backgrounds due to its large production rate, especially
when the value of $\mu$ is small. Another type of dominant backgrounds are the
QCD multi-jets plus $W$ or $Z$ boson. In the case of multi-jets plus a $W$
boson, it can fake the signal when the $W$ boson decays to $\ell\nu$ with the
charged lepton missing detection (too soft, travelling along the beamline or
too close to a jet) or the $W$ boson decaying to $\tau\nu$ with the secondary
jet from the hadronic $\tau$ decay missing detection. For the case of multi-
jets plus a $Z$ boson, it can fake the signal when $Z$ decays to
$\nu\bar{\nu}$. Although their cross sections are about 1-2 pb, mucher larger
than the cross section of the signal, they are not the major backgrounds due
to the small fake efficiency. The $t\bar{t}h(h\to b\bar{b}),t\bar{t}b\bar{b}$
process, which were not the main background in the preceding section, is now
likely the major background. So we consider the backgrounds from $bbbb$,
$bbbbZ$, $bbjjW^{\pm}(W^{\pm}\rightarrow\ell\nu)$,
$bbjjW^{\pm}(W^{\pm}\rightarrow\tau\nu)$, $bbjjZ$ and $t\bar{t}h(h\to
b\bar{b}),t\bar{t}b\bar{b}$ processes. All these backgrounds are generated
using MadGraph with the PDF CTEQ6L1 2003JPP  at the leading order. Since the
cross-section of this signal is large enough even for a relatively large
$\mu$, we will choose different values of $\mu$ to show the results. At the
parton level, we simply require $p_{Tj}>30$ GeV, $|\eta_{j}|<2.7$, $\Delta
R_{jj}>0.4$ in our simulations.

We note that the triggers and jet reconstructions are very important for the
signal, in which the transverse momentum of the Higgs boson plays a key role.
So after parton-shower and hadronization, we show the $P_{T}$ of the leading
Higgs and ${\>/\\!\\!\\!\\!E_{T}}$ distributions in Fig. 4. We find that the
signal has relatively low transverse momentum for the Higgs, which is
typically below 200 GeV for the value of $\mu$ below 400 GeV. We require at
least four anti-$k_{t}$ R=0.4 $b-$jets for the signal (the efficiency for
reconstructing a Higgs boson from two anti-$k_{t}$ $R=0.4$ jets is higher than
reconstructing from a Cambridge-Aachen jet with $R=1.2$ in this case 2015NHP
).

Figure 4:  Distributions of $P_{T}$ of the leading Higgs and
${\>/\\!\\!\\!\\!E_{T}}$ for the signal
$bb\bar{b}\bar{b}+{\>/\\!\\!\\!\\!E_{T}}$ with different values of $\mu$. The
$y$-axis represents the normalized number of events.

So the four $b-$tagged jets are required, paired into two dijets, to
reconstruct the Higgs bosons, and this is a powerful way to reduce the
backgrounds.Obviously, there are three different ways for each event to pair
the b jets. Here the traditional $\chi^{2}$ method is used to pick out the
reconstruction which has the minimum value of $\chi^{2}$.

Figure 5:  Distributions of signal and backgrounts after basic cuts for
$\mu=150$ GeV. The $y$-axis represents the normalized number of events

The distributions of the signal and backgrounds are displayed in Fig. 5 for
$\mu=150$ GeV. It is easy to find that the distributions of the invariant
masses of the leading and next-to-leading Higgs bosons ($M_{HL}$ and $M_{HN}$)
and the corresponding azimuths ($\Delta R_{HL}$ and $\Delta R_{HN}$) could be
used to effectively suppress the background. The cuts 120 GeV $<M_{H(N,L)}<$
130 GeV and $\Delta R_{H(N,L)}<2$ are used for the two reconstructed Higgs
bosons. Note that the rather strong invariant mass cuts are applied since the
$bbbb$ background is large. Here the missing energy ${\>/\\!\\!\\!\\!E_{T}}$
cannot suppress the QCD background very well, and thus we impose a relatively
strong cut on missing energy ${\>/\\!\\!\\!\\!E_{T}}>$ 50 GeV to get a larger
value of $S/B$. With all above cuts, the significance can only reach
$0.679\sigma$, as shown in Table 2.

Table 2:  Signal and background cross sections in fb after the cuts for $\mu$
= 150 GeV. The significance $\sigma$ is calculated for a luminosity of 3
$\rm{ab}^{-1}$ at the 14 TeV LHC.

Cuts | signal | $bbbbZ$ | $bbjjW_{1}$ | $bbjjW_{2}$ | $bbjjZ$ | $tth(b\bar{b})$ | $b\bar{b}b\bar{b}$ | $t\bar{t}b\bar{b}$ | $\sigma$ | $S/B$
---|---|---|---|---|---|---|---|---|---|---
130 GeV > $M_{H(L,N)}$ > 120 GeV | 0.036 | 0.0014 | 0.005 | 0.003 | 0.005 | 0.057 | 22.41 | 0.086 | 0.53 | 0.0016
$\Delta R_{H(L,N)}$ < 2 | 0.017 | 0 | 0.001 | 0.0005 | 0 | 0.028 | 3.58 | 0.04 | 0.637 | 0.005
${\>/\\!\\!\\!\\!E_{T}}$ > 50 GeV | 0.0073 | 0 | 0.001 | 0 | 0 | 0.028 | 0.54 | 0.02 | 0.679 | 0.012

With the increase of $\mu$, the cross section of the signal decreases rapidly,
and at the same time the missing energy and $P_{T}^{H}$ become larger, as
shown in Fig. 4. Now we take $\mu$ = 300 GeV as an example to show the
analysis for a larger $\mu$. The cuts 105 GeV <$M_{H(L,N)}$ <140 GeV, which
keep most of the signal events, are used for paired b jets. Then, since there
are no obvious differences on angular separations between signal and
backgrounds, we use the same cuts $\Delta R_{H(N,L)}<2$ as before. Unlike the
case of $\mu$ = 150 GeV, in this case the missing energy
${\>/\\!\\!\\!\\!E_{T}}$ could play an important role in distinguishing signal
and backgrounds. We set ${\>/\\!\\!\\!\\!E_{T}}>190$ GeV and this cut can
increase the significance by about a factor of ten. The detailed results could
be found in Table 3.

Table 3:  Signal and background cross sections in fb after cuts at different
stages for $\mu$ = 300 GeV.

Cuts | signal | $bbbbZ$ | $bbjjW_{1}$ | $bbjjW_{2}$ | $bbjjZ$ | $tth(b\bar{b})$ | $b\bar{b}b\bar{b}$ | $t\bar{t}b\bar{b}$ | $\sigma$ | $S/B$
---|---|---|---|---|---|---|---|---|---|---
140 GeV > $M_{H(L,N)}$ > 105 GeV | 0.039 | 0.018 | 0.044 | 0.024 | 0.054 | 0.6 | 236.442 | 0.75 | 0.138 | 0.00013
$\Delta R_{H(L,N)}$ < 2 | 0.028 | 0.0047 | 0.014 | 0.007 | 0.014 | 0.31 | 47.53 | 0.39 | 0.21 | 0.00058
${\>/\\!\\!\\!\\!E_{T}}$ > 190 GeV | 0.0079 | 0.0005 | 0.001 | 0.0005 | 0.0027 | 0.011 | 0.0173 | 0.011 | 2.02 | 0.18

From the above results we find that the ${\>/\\!\\!\\!\\!E_{T}}$ cut plays a
major role in distinguishing the signal from the backgrounds for a large
$\mu$. The detailed cuts are summarized as follows:

  * •

The four leading jets must be $b-$tagged with $P_{T}$ > 40 GeV;

  * •

Exactly two Higgs are reconstructed, with 140 GeV >$M_{H(L,N)}$ > 105 GeV;

  * •

Proximity cut of $\Delta R_{H(N,L)}$ < 2 for the two Higgs;

  * •

${\>/\\!\\!\\!\\!E_{T}}>190$ GeV for the reconstructed missing transverse
momentum.

After the sequential cuts, we find that the multi-jets background are greatly
suppressed, while the $t\bar{t}h,b\bar{b}b\bar{b},t\bar{t}b\bar{b}$ background
remains large and will be further cut with a larger ${\>/\\!\\!\\!\\!E_{T}}$.
With the increase of $\mu$ value, the signal cross-section will be reduced but
the cut with a larger ${\>/\\!\\!\\!\\!E_{T}}$ can help for the signal
significance. For $\mu=300$ GeV, the signal significance can reach to
2.02$\sigma$, as shown in Fig. 6. As $\mu$ increases above 320 GeV the
significance drops below 2$\sigma$. Note that in this scenario the $\mu$ value
determines the lightest neutralino mass.

Figure 6:  The signal $bb\bar{b}\bar{b}+{\>/\\!\\!\\!\\!E_{T}}$ significance
versus the $\mu$ value at the HL-LHC (3 ab-1, 14 TeV).

###  3.3 The signal of $hhG^{\prime}G^{\prime}\rightarrow
b\bar{b}\gamma\gamma+{\>/\\!\\!\\!\\!E_{T}}$

Finally we turn to the third channel with $h\rightarrow b\bar{b}$ and
$h\rightarrow\gamma\gamma$. It has a pair of photons reconstructed at the
invariant mass around the Higgs boson. So the backgrounds include the single
Higgs production such as $t\bar{t}H$. For non-resonant backgrounds and jet-
fake backgrounds, we should consider missing energy based on the backgrounds
considered in the preceding section. From the calculation we mainly need to
consider the backgounds from $bbj\gamma
W^{\pm}(W^{\pm}\rightarrow\ell^{\pm}\nu)$, $bbj\gamma
W^{\pm}(W^{\pm}\rightarrow\tau^{\pm}\nu)$, $bbj\gamma
Z(Z\rightarrow\nu\bar{\nu})$, $bbbbZ$, and $t\bar{t}\gamma$ productions.

The parton-level events for backgrounds are generated with MadGraph5, using
the PDF CTEQ6L1 2003JPP . For the signal calculation, it involves the loop-
induced Higgs decay to di-photon and we use the loop-induced model
Hirschi:2015iia . The parton-level cuts are imposed in order to avoid any
divergence in the parton-level calculation ATLAS:2017muo : $P_{Tj}>20$ GeV,
$P_{T\gamma}>25$ GeV, $|\eta_{j}|<2.5$, $|\eta_{\gamma}|<2.5$, $m_{jj}>25$
GeV.

Figure 7:  Distributions for the signal with $\mu$ =150 GeV and backgrounts.
The $y$-axis represents the normalized number of events for each process.

For signal and backgrounds we require a pair of isolated photons and a pair of
isolated $b$-tagged jets which are reconstructed near the Higgs boson mass. In
Fig. 7, we show the distributions of $P_{T}^{bb}$, $P_{T}^{\gamma\gamma}$,
${\>/\\!\\!\\!\\!E_{T}}$, $m_{\gamma\gamma}$ and $m_{bb}$ for signal and
backgrounds. Since the angular separations have no obvious difference between
signal and backgrounds, we set loose cuts $1.0<\Delta R_{bb}<3.8$ and
$1<\Delta R_{\gamma\gamma}<3.7$. Obviously, $m_{\gamma\gamma}$ and
$P_{T}^{\gamma\gamma}$ have distinctive features to distinguish signal from
backgrounds expect $t\bar{t}h$. In order to reduce the background of
$t\bar{t}h$, we require $P_{T}^{bb}$ in the range of $[80,250]$ GeV and
$P_{T}^{\gamma\gamma}>70$ GeV. The ${\>/\\!\\!\\!\\!E_{T}}$ cut can
efficiently suppress the multi-jet backgrounds and we set 10 GeV
$<{\>/\\!\\!\\!\\!E_{T}}<90$ GeV. The cuts can be summarized in the following:

  * •

The two leading jets must be $b-$tagged with 250 GeV $>p_{T}^{bb}>80$ GeV;

  * •

Exactly two isolated photons with $p_{T}^{\gamma\gamma}>70$ GeV;

  * •

Proximity cut of $1<\Delta R{{}_{\gamma\gamma}}$<3.7 for the two photons;

  * •

Proximaty cut of $1<\Delta R{{}_{bb}}<3.8$ for the two $b-$tagged jets;

  * •

121 GeV $<m_{\gamma\gamma}<128$ GeV for the two photons;

  * •

100 GeV $<m_{bb}<150$ GeV for the two $b-$tagged jets;

  * •

90 GeV $>{\>/\\!\\!\\!\\!E_{T}}=|{\;/\\!\\!\\!\\!\vec{P}_{T}}|>25$ GeV for the
reconstructed missing transverse momentum.

The detailed results are shown in Table 4. From our analysis we find that the
$t\bar{t}\gamma$ background still remains large after these cuts. The signal
significance can only reach to $0.79\sigma$.

Table 4:  The cut flow for signal and background cross sections in fb for a
luminosity of 3 $\rm{ab}^{-1}$ at the 14 TeV LHC.

Cuts | signal | $b\bar{b}j\gamma W_{1}$ | $b\bar{b}j\gamma W_{2}$ | $b\bar{b}jZ$ | $t\bar{t}\gamma$ | $t\bar{t}h$ | $\sigma$ | $S/B$
---|---|---|---|---|---|---|---|---
80 GeV < $P_{T}^{bb}$ < 250 GeV, $P_{T}^{\gamma\gamma}$ >70 GeV | 0.0093 | 0.0466 | 0.036 | 0.044 | 4.09 | 0.011 | 0.25 | 0.0022
100 GeV < $m_{bb}$ < 150 GeV, GeV 121 $<m_{\gamma\gamma}$ < 128 GeV | 0.0051 | 0.0008 | 0.001 | 0.0005 | 0.089 | 0.0026 | 0.9 | 0.054
90 GeV > ${\>/\\!\\!\\!\\!E_{T}}$ > 25 GeV | 0.0035 | 0.0006 | 0.0009 | 0 | 0.052 | 0.002 | 0.79 | 0.62

Some remarks are in order. From the above results we know that for the three
channels only the $bbbb+{\>/\\!\\!\\!\\!E_{T}}$ signal might be observable at
the HL-LHC. Note that some machine learning methods Blanke:2019hpe ;
Amacker:2020bmn ; Abdughani:2019wuv ; Abdughani:2020xfo  which have been used
for many body final states, can help to enhance the signal significance. Since
our main aim is to show the new signature that comes up with pseudo-goldstino
and to study the role of ${\>/\\!\\!\\!\\!E_{T}}$ in distinguishing the signal
from backgound, only the traditional cut-flow Monte Carlo simulations were
used in our analysis. In addition, we did not consider the explanation of the
muon $g-2$ in our scenario, which was found (see, e.g., Abdughani:2019wai ;
Athron:2021iuf ; Wang:2021bcx ; Endo:2021zal ) to require light
neutralinos/charginos and light sleptons, possibly accessible at the HL-LHC
Abdughani:2019wai . For the status of low-energy supersymmetry confronted with
current experimental constraints, such as collider searches, dark matter
searches and muon $g-2$ measurement, recent overviews could be found in Refs.
Baer:2020kwz ; Wang:2022rfd .

Finally, we point out that our study may be also applicable to higgsino decay
into a light singlino (which can be rather light as the lightest superparticle
Cao:2013mqa ) plus a Higgs boson in the NMSSM Ellwanger:2009dp  which
introduces a singlet Higgs superfield, mixing with the Higgs doublets, and
thus relieves the 125 GeV Higgs mass constraint on stop masses Cao:2012fz .

##  4 Conclusion

We considered the scenario of multi-sector SUSY breaking which predicts light
pseudo-goldstinos and opens the decay mode of neutralino into pseudo-goldstino
plus Higgs boson insider the detector at the LHC, leading to the signal of
Higgs pair plus missing energy. We studied the observability of such Higgs
pair plus missing energy from the decay of neutralino produced at the HL-LHC
(14 TeV with a luminosity of 3 $\rm{ab}^{-1}$). We considered light higgsinos
assumed in natural SUSY as an example and studied the three decay channels of
the Higgs pair ($bbWW^{*}$, $bb\gamma\gamma$, $bbbb$). From detailed Monte
Carlo simulations for the signal and backgrounds, the following observations
are obtained: (i) The best channel is $bbbb+{\>/\\!\\!\\!\\!E_{T}}$, whose
statistical significance can reach $2\sigma$ for a light higgsino in natural
SUSY allowed by current experiments; (ii) The channels
$bb\gamma\gamma+{\>/\\!\\!\\!\\!E_{T}}$ and $bbWW^{*}+{\>/\\!\\!\\!\\!E_{T}}$
respectively give maximal statistical significance of $1\sigma$ and
$0.79\sigma$. So the Higgs pair plus missing energy signal from the higgsino
decay can be over the SM Higgs pair result which is about $1.8\sigma$ at the
HL-LHC.

## Acknowledgments

This work was supported in part by IHEP under Grant No. Y9515570U1, by the
National Natural Science Foundation of China (NNSFC) under grant Nos. 11821505
and 12075300, by Peng-Huan-Wu Theoretical Physics Innovation Center
(12047503), by the CAS Center for Excellence in Particle Physics (CCEPP), by
the CAS Key Research Program of Frontier Sciences, and by a Key R&D Program of
Ministry of Science and Technology of China under number 2017YFA0402204, and
by the Key Research Program of the Chinese Academy of Sciences, Grant NO.
XDPB15.

## Note added:

When preparing this manuscript, the CMS collaboration reported their search
results for the signal of Higgs pair plus missing energy using an integrated
luminosity of 137 fb-1 CMS:2022vpy . When applied to higgsino decay into Higgs
plus pseudo-goldstino, a range 305 GeV $>m_{\tilde{\chi}_{1}^{0}}(\simeq\mu)>$
265 GeV can be excluded by this CMS search, assuming the decay branching ratio
to be 1. Since in our study we assumed the branching ratio of higgsino decay
to Higgs plus pseudo-goldstino to be 0.5, this excluded range 305 GeV
$>m_{\tilde{\chi}_{1}^{0}}(\simeq\mu)>$ 265 GeV will become much narrower when
applied to our scenario.

## References

  * (1) ATLAS collaboration, M. Aaboud et al., _Search for Higgs boson pair production in the $b\bar{b}WW^{*}$ decay mode at $\sqrt{s}=13$ TeV with the ATLAS detector_, 1811.04671.
  * (2) CMS collaboration, _Search for resonant and non-resonant Higgs boson pair production in the $\mathrm{b}\overline{\mathrm{b}}\mathrm{l}\nu\mathrm{l}\nu$ final state at $\sqrt{s}=13~{}\mathrm{TeV}$_, CMS-PAS-HIG-17-006.
  * (3) ATLAS collaboration, M. Aaboud et al., _Search for pair production of Higgs bosons in the $b\bar{b}b\bar{b}$ final state using proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _JHEP_ 01 (2019) 030, [1804.06174].
  * (4) CMS collaboration, _Search for Non-Resonant Higgs Pair-Production in the $\mathrm{b}\overline{\mathrm{b}}\mathrm{b}\overline{\mathrm{b}}$ Final State with the CMS detector_, CMS-PAS-HIG-17-017.
  * (5) CMS collaboration, A. M. Sirunyan et al., _Search for Higgs boson pair production in the $\gamma\gamma\mathrm{b\overline{b}}$ final state in pp collisions at $\sqrt{s}=$ 13 TeV_, _Phys. Lett._ B788 (2019) 7–36, [1806.00408].
  * (6) ATLAS collaboration, M. Aaboud et al., _Search for Higgs boson pair production in the $\gamma\gamma b\bar{b}$ final state with 13 TeV $pp$ collision data collected by the ATLAS experiment_, _JHEP_ 11 (2018) 040, [1807.04873].
  * (7) ATLAS collaboration, M. Aaboud et al., _Search for resonant and non-resonant Higgs boson pair production in the ${b\bar{b}\tau^{+}\tau^{-}}$ decay channel in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Phys. Rev. Lett._ 121 (2018) 191801, [1808.00336].
  * (8) CMS collaboration, A. M. Sirunyan et al., _Search for Higgs boson pair production in events with two bottom quarks and two tau leptons in proton-proton collisions at $\sqrt{s}$ =13 TeV_, _Phys. Lett._ B778 (2018) 101–127, [1707.02909].
  * (9) ATLAS collaboration, M. Aaboud et al., _Search for Higgs boson pair production in the $WW^{(*)}WW^{(*)}$ decay channel using ATLAS data recorded at $\sqrt{s}=13$ TeV_, 1811.11028.
  * (10) G. Aad, B. Abbott, D. Abbott, A. Abed Abud, K. Abeling, D. Abhayasinghe et al., _Combination of searches for higgs boson pairs in pp collisions at $\sqrt{s}$ = 13 tev with the atlas detector_, _Physics Letters B_ 800 (Jan, 2020) 135103.
  * (11) R. Contino, C. Grojean, M. Moretti, F. Piccinini and R. Rattazzi, _Strong Double Higgs Production at the LHC_ , _JHEP_ 05 (2010) 089, [1002.1011].
  * (12) R. Grober and M. Muhlleitner, _Composite Higgs Boson Pair Production at the LHC_ , _JHEP_ 06 (2011) 020, [1012.1562].
  * (13) M. J. Dolan, C. Englert and M. Spannowsky, _New Physics in LHC Higgs boson pair production_ , _Phys. Rev. D_ 87 (2013) 055002, [1210.8166].
  * (14) J. Cao, Z. Heng, L. Shang, P. Wan and J. M. Yang, _Pair Production of a 125 GeV Higgs Boson in MSSM and NMSSM at the LHC_ , _JHEP_ 04 (2013) 134, [1301.6437].
  * (15) C. Han, X. Ji, L. Wu, P. Wu and J. M. Yang, _Higgs pair production with SUSY QCD correction: revisited under current experimental constraints_ , _JHEP_ 04 (2014) 003, [1307.3790].
  * (16) B. Hespel, D. Lopez-Val and E. Vryonidou, _Higgs pair production via gluon fusion in the Two-Higgs-Doublet Model_ , _JHEP_ 09 (2014) 124, [1407.0281].
  * (17) J. Cao, D. Li, L. Shang, P. Wu and Y. Zhang, _Exploring the Higgs Sector of a Most Natural NMSSM and its Prediction on Higgs Pair Production at the LHC_ , _JHEP_ 12 (2014) 026, [1409.8431].
  * (18) S. Dawson, A. Ismail and I. Low, _What’s in the loop? The anatomy of double Higgs production_ , _Phys. Rev. D_ 91 (2015) 115008, [1504.05596].
  * (19) L.-C. Lü, C. Du, Y. Fang, H.-J. He and H. Zhang, _Searching heavier Higgs boson via di-Higgs production at LHC Run-2_ , _Phys. Lett. B_ 755 (2016) 509–522, [1507.02644].
  * (20) Z. Kang, P. Ko and J. Li, _New Physics Opportunities in the Boosted Di-Higgs-Boson Plus Missing Transverse Energy Signature_ , _Phys. Rev. Lett._ 116 (2016) 131801, [1504.04128].
  * (21) Q.-H. Cao, B. Yan, D.-M. Zhang and H. Zhang, _Resolving the degeneracy in single higgs production with higgs pair production_ , _Physics Letters B_ 752 (Jan, 2016) 285–290.
  * (22) Q.-H. Cao, Y. Liu and B. Yan, _Measuring trilinear higgs coupling in whh and zhh productions at the high-luminosity lhc_ , _Physical Review D_ 95 (Apr, 2017) .
  * (23) Q.-H. Cao, G. Li, B. Yan, D.-M. Zhang and H. Zhang, _Double higgs production at the 14 tev lhc and a 100 tev pp collider_ , _Physical Review D_ 96 (Nov, 2017) .
  * (24) K. Nakamura, K. Nishiwaki, K.-y. Oda, S. C. Park and Y. Yamamoto, _Di-higgs enhancement by neutral scalar as probe of new colored sector_ , _Eur. Phys. J. C_ 77 (2017) 273, [1701.06137].
  * (25) J. Chang, C.-R. Chen and C.-W. Chiang, _Higgs boson pair productions in the Georgi-Machacek model at the LHC_ , _JHEP_ 03 (2017) 137, [1701.06291].
  * (26) J. Ren, R.-Q. Xiao, M. Zhou, Y. Fang, H.-J. He and W. Yao, _LHC Search of New Higgs Boson via Resonant Di-Higgs Production with Decays into 4W_ , _JHEP_ 06 (2018) 090, [1706.05980].
  * (27) P. Huang, A. Joglekar, M. Li and C. E. M. Wagner, _Corrections to di-Higgs boson production with light stops and modified Higgs couplings_ , _Phys. Rev. D_ 97 (2018) 075001, [1711.05743].
  * (28) A. Adhikary, S. Banerjee, R. K. Barman, B. Bhattacherjee and S. Niyogi, _Revisiting the non-resonant Higgs pair production at the HL-LHC_ , _JHEP_ 07 (2018) 116, [1712.05346].
  * (29) P. Agrawal, D. Saha, L.-X. Xu, J.-H. Yu and C. P. Yuan, _Determining the shape of the Higgs potential at future colliders_ , _Phys. Rev. D_ 101 (2020) 075023, [1907.02078].
  * (30) K. Cheung, A. Jueid, C.-T. Lu, J. Song and Y. W. Yoon, _Disentangling new physics effects on nonresonant Higgs boson pair production from gluon fusion_ , _Phys. Rev. D_ 103 (2021) 015019, [2003.11043].
  * (31) G. Li, L.-X. Xu, B. Yan and C.-P. Yuan, _Resolving the degeneracy in top quark yukawa coupling with higgs pair production_ , _Physics Letters B_ 800 (Jan, 2020) 135070.
  * (32) H. Abouabid, A. Arhrib, D. Azevedo, J. E. Falaki, P. M. Ferreira, M. Mühlleitner et al., _Benchmarking Di-Higgs Production in Various Extended Higgs Sector Models_ , 2112.12515.
  * (33) L. Wu, J. M. Yang, C.-P. Yuan and M. Zhang, _Higgs self-coupling in the MSSM and NMSSM after the LHC Run 1_ , _Phys. Lett. B_ 747 (2015) 378–389, [1504.06932].
  * (34) L. Wang, W. Wang, J. M. Yang and H. Zhang, _Higgs-pair production in littlest Higgs model with T-parity_ , _Phys. Rev. D_ 76 (2007) 017702, [0705.3392].
  * (35) X.-F. Han, L. Wang and J. M. Yang, _Higgs-pair Production and Decay in Simplest Little Higgs Model_ , _Nucl. Phys. B_ 825 (2010) 222–230, [0908.1827].
  * (36) R. Argurio, Z. Komargodski and A. Mariotti, _Pseudo-Goldstini in Field Theory_ , _Phys. Rev. Lett._ 107 (2011) 061601, [1102.2386].
  * (37) J. Dai, T. Liu and J. M. Yang, _An explicit calculation of pseudo-goldstino mass at the leading three-loop level_ , _JHEP_ 06 (2021) 175, [2104.12656].
  * (38) C. Cheung, Y. Nomura and J. Thaler, _Goldstini_ , _JHEP_ 03 (2010) 073, [1002.1967].
  * (39) C. Cheung, J. Mardon, Y. Nomura and J. Thaler, _A Definitive Signal of Multiple Supersymmetry Breaking_ , _JHEP_ 07 (2010) 035, [1004.4637].
  * (40) N. Craig, J. March-Russell and M. McCullough, _The Goldstini Variations_ , _JHEP_ 10 (2010) 095, [1007.1239].
  * (41) J. Thaler and Z. Thomas, _Goldstini Can Give the Higgs a Boost_ , _JHEP_ 07 (2011) 060, [1103.1631].
  * (42) C. Cheung, F. D’Eramo and J. Thaler, _The Spectrum of Goldstini and Modulini_ , _JHEP_ 08 (2011) 115, [1104.2600].
  * (43) R. Argurio, K. De Causmaecker, G. Ferretti, A. Mariotti, K. Mawatari and Y. Takaesu, _Collider signatures of goldstini in gauge mediation_ , _JHEP_ 06 (2012) 096, [1112.5058].
  * (44) E. Dudas, G. von Gersdorff, D. M. Ghilencea, S. Lavignac and J. Parmentier, _On non-universal Goldstino couplings to matter_ , _Nucl. Phys. B_ 855 (2012) 570–591, [1106.5792].
  * (45) T. Liu, L. Wang and J. M. Yang, _Higgs decay to goldstini and its observability at the LHC_ , _Phys. Lett. B_ 726 (2013) 228–233, [1301.5479].
  * (46) K.-i. Hikasa, T. Liu, L. Wang and J. M. Yang, _Pseudo-goldstino and electroweak gauginos at the LHC_ , _JHEP_ 07 (2014) 065, [1403.5731].
  * (47) T. Liu, L. Wang and J. M. Yang, _Pseudo-goldstino and electroweakinos via VBF processes at LHC_ , _JHEP_ 02 (2015) 177, [1411.6105].
  * (48) G. Ferretti, A. Mariotti, K. Mawatari and C. Petersson, _Multiphoton signatures of goldstini at the LHC_ , _JHEP_ 04 (2014) 126, [1312.1698].
  * (49) D. B. Franzosi, G. Ferretti, E. Riefel and S. Strandberg, _Electroweak signatures of gauge-mediated supersymmetry breaking in multiple hidden sectors_ , _JHEP_ 01 (2022) 139, [2111.04775].
  * (50) C. Brust, A. Katz, S. Lawrence and R. Sundrum, _SUSY, the Third Generation and the LHC_ , _JHEP_ 03 (2012) 103, [1110.6670].
  * (51) M. Papucci, J. T. Ruderman and A. Weiler, _Natural SUSY Endures_ , _JHEP_ 09 (2012) 035, [1110.6926].
  * (52) L. J. Hall, D. Pinner and J. T. Ruderman, _A Natural SUSY Higgs Near 126 GeV_ , _JHEP_ 04 (2012) 131, [1112.2703].
  * (53) J. L. Feng and D. Sanford, _A Natural 125 GeV Higgs Boson in the MSSM from Focus Point Supersymmetry with A-Terms_ , _Phys. Rev. D_ 86 (2012) 055015, [1205.2372].
  * (54) H. Baer, V. Barger, P. Huang, A. Mustafayev and X. Tata, _Radiative natural SUSY with a 125 GeV Higgs boson_ , _Phys. Rev. Lett._ 109 (2012) 161802, [1207.3343].
  * (55) X. Tata, _Natural supersymmetry: status and prospects_ , _Eur. Phys. J. ST_ 229 (2020) 3061–3083, [2002.04429].
  * (56) C. Han, A. Kobakhidze, N. Liu, A. Saavedra, L. Wu and J. M. Yang, _Probing Light Higgsinos in Natural SUSY from Monojet Signals at the LHC_ , _JHEP_ 02 (2014) 049, [1310.4274].
  * (57) GAMBIT collaboration, P. Athron et al., _Combined collider constraints on neutralinos and charginos_ , _Eur. Phys. J. C_ 79 (2019) 395, [1809.02097].
  * (58) B. Allanach, _Softsusy: A program for calculating supersymmetric spectra_ , _Computer Physics Communications_ 143 (Mar, 2002) 305–331.
  * (59) S. Borowka, N. Greiner, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk et al., _Higgs Boson Pair Production in Gluon Fusion at Next-to-Leading Order with Full Top-Quark Mass Dependence_ , _Phys. Rev. Lett._ 117 (2016) 012001, [1604.06447].
  * (60) J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, M. Spira and J. Streicher, _Gluon fusion into Higgs pairs at NLO QCD and the top mass scheme_ , _Eur. Phys. J. C_ 79 (2019) 459, [1811.05692].
  * (61) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ , _JHEP_ 07 (2014) 079, [1405.0301].
  * (62) NNPDF collaboration, R. D. Ball, V. Bertone, S. Carrazza, L. Del Debbio, S. Forte, A. Guffanti et al., _Parton distributions with QED corrections_ , _Nucl. Phys._ B877 (2013) 290–320, [1308.0598].
  * (63) W. Beenakker, M. Klasen, M. Kramer, T. Plehn, M. Spira and P. M. Zerwas, _The Production of charginos / neutralinos and sleptons at hadron colliders_ , _Phys. Rev. Lett._ 83 (1999) 3780–3783, [hep-ph/9906298].
  * (64) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, _Feynrules 2.0 — a complete toolbox for tree-level phenomenology_ , _Computer Physics Communications_ 185 (Aug, 2014) 2250–2300.
  * (65) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer and T. Reiter, _Ufo – the universal feynrules output_ , _Computer Physics Communications_ 183 (Jun, 2012) 1201–1214.
  * (66) T. Sjostrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An Introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159–177, [1410.3012].
  * (67) DELPHES 3 collaboration, J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaitre, A. Mertens et al., _DELPHES 3, A modular framework for fast simulation of a generic collider experiment_ , _JHEP_ 02 (2014) 057, [1307.6346].
  * (68) M. Cacciari, G. P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J._ C72 (2012) 1896, [1111.6097].
  * (69) M. Cacciari, G. P. Salam and G. Soyez, _The anti- $k_{t}$ jet clustering algorithm_, _JHEP_ 04 (2008) 063, [0802.1189].
  * (70) ATLAS collaboration, _Expected performance of the ATLAS detector at the High-Luminosity LHC_ , Tech. Rep. ATL-PHYS-PUB-2019-005, CERN, Geneva, Jan, 2019.
  * (71) J. H. Kim, M. Kim, K. Kong, K. T. Matchev and M. Park, _Portraying double higgs at the large hadron collider_ , _Journal of High Energy Physics_ 2019 (Sep, 2019) .
  * (72) M. Czakon, P. Fiedler and A. Mitov, _Total Top-Quark Pair-Production Cross Section at Hadron Colliders Through ${\cal O}(\alpha_{S}^{4})$_, _Phys. Rev. Lett._ 110 (2013) 252004, [1303.6254].
  * (73) LHC Higgs Cross Section Working Group collaboration, S. Dittmaier et al., _Handbook of LHC Higgs Cross Sections: 1. Inclusive Observables_ , 1101.0593.
  * (74) LHC Higgs Cross Section Working Group collaboration, D. de Florian et al., _Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector_ , 1610.07922.
  * (75) D. Wardrope, E. Jansen, N. Konstantinidis, B. Cooper, R. Falla and N. Norjoharuddeen, _Non-resonant higgs-pair production in the $b\bar{b}b\bar{b}$ final state at the lhc_, _The European Physical Journal C_ 75 (May, 2015) .
  * (76) D. Stump, J. Huston, J. Pumplin, W.-K. Tung, H.-L. Lai, S. Kuhlmann et al., _Inclusive jet production, parton distributions, and the search for new physics_ , _Journal of High Energy Physics_ 2003 (Oct, 2003) 046–046.
  * (77) V. Hirschi and O. Mattelaer, _Automated event generation for loop-induced processes_ , _JHEP_ 10 (2015) 146, [1507.00020].
  * (78) ATLAS collaboration, _Study of the double Higgs production channel $H(\rightarrow b\bar{b})H(\rightarrow\gamma\gamma)$ with the ATLAS experiment at the HL-LHC_, .
  * (79) M. Blanke, S. Kast, J. M. Thompson, S. Westhoff and J. Zurita, _Spotting hidden sectors with Higgs binoculars_ , _JHEP_ 04 (2019) 160, [1901.07558].
  * (80) J. Amacker et al., _Higgs self-coupling measurements using deep learning in the $b\overline{b}b\overline{b}$ final state_, _JHEP_ 12 (2020) 115, [2004.04240].
  * (81) M. Abdughani, J. Ren, L. Wu, J. M. Yang and J. Zhao, _Supervised deep learning in high energy phenomenology: a mini review_ , _Commun. Theor. Phys._ 71 (2019) 955, [1905.06047].
  * (82) M. Abdughani, D. Wang, L. Wu, J. M. Yang and J. Zhao, _Probing the triple Higgs boson coupling with machine learning at the LHC_ , _Phys. Rev. D_ 104 (2021) 056003, [2005.11086].
  * (83) M. Abdughani, K.-I. Hikasa, L. Wu, J. M. Yang and J. Zhao, _Testing electroweak SUSY for muon $g$ $-$ 2 and dark matter at the LHC and beyond_, _JHEP_ 11 (2019) 095, [1909.07792].
  * (84) P. Athron, C. Balázs, D. H. J. Jacob, W. Kotlarski, D. Stöckinger and H. Stöckinger-Kim, _New physics explanations of a μ in light of the FNAL muon g-2 measurement_, _JHEP_ 09 (2021) 080, [2104.03691].
  * (85) F. Wang, L. Wu, Y. Xiao, J. M. Yang and Y. Zhang, _GUT-scale constrained SUSY in light of new muon g-2 measurement_ , _Nucl. Phys. B_ 970 (2021) 115486, [2104.03262].
  * (86) M. Endo, K. Hamaguchi, S. Iwamoto and T. Kitahara, _Supersymmetric interpretation of the muon g – 2 anomaly_ , _JHEP_ 07 (2021) 075, [2104.03217].
  * (87) H. Baer, V. Barger, S. Salam, D. Sengupta and K. Sinha, _Status of weak scale supersymmetry after LHC Run 2 and ton-scale noble liquid WIMP searches_ , _Eur. Phys. J. ST_ 229 (2020) 3085–3141, [2002.03013].
  * (88) F. Wang, W. Wang, J. M. Yang, Y. Zhang and B. Zhu, _Low energy supersymmetry confronted with current experiments: an overview_ , 2201.00156.
  * (89) J. Cao, C. Han, L. Wu, P. Wu and J. M. Yang, _A light SUSY dark matter after CDMS-II, LUX and LHC Higgs data_ , _JHEP_ 05 (2014) 056, [1311.0678].
  * (90) U. Ellwanger, C. Hugonie and A. M. Teixeira, _The Next-to-Minimal Supersymmetric Standard Model_ , _Phys. Rept._ 496 (2010) 1–77, [0910.1785].
  * (91) J.-J. Cao, Z.-X. Heng, J. M. Yang, Y.-M. Zhang and J.-Y. Zhu, _A SM-like Higgs near 125 GeV in low energy SUSY: a comparative study for MSSM and NMSSM_ , _JHEP_ 03 (2012) 086, [1202.5821].
  * (92) CMS collaboration, A. Tumasyan et al., _Search for higgsinos decaying to two Higgs bosons and missing transverse momentum in proton-proton collisions at $\sqrt{s}$ = 13 TeV_, 2201.04206.