# Relationship between Successive Flares in the Same Active Region and Space-
Weather HMI Active Region Patch (SHARP) Parameters

Hao Ran  State Key Laboratory of Space Weather, National Space Science Center,
Chinese Academy of Sciences, Beijing 100190, China<EMAIL_ADDRESS>University of Chinese Academy of Sciences, Beijing 100049, China  Ying D. Liu
State Key Laboratory of Space Weather, National Space Science Center, Chinese
Academy of Sciences, Beijing 100190, China<EMAIL_ADDRESS>University of
Chinese Academy of Sciences, Beijing 100049, China  Yang Guo  School of
Astronomy and Space Science, Nanjing University, Nanjing 210023, China  Rui
Wang  State Key Laboratory of Space Weather, National Space Science Center,
Chinese Academy of Sciences, Beijing 100190, China<EMAIL_ADDRESS>
###### Abstract

A solar active region (AR) may produce multiple notable flares during its
passage across the solar disk. We investigate successive flares from flare-
eruptive active regions, and explore their relationship with solar magnetic
parameters. We examine six ARs in this study, each with at least one major
flare above X1.0. The Space-Weather HMI Active Region Patch (SHARP) is
employed in this study to parameterize the ARs. We aim to identify the most
flare-related SHARP parameters and lay foundation for future practical flare
forecasts. We first evaluate the correlation coefficients between the SHARP
parameters and the successive flare production. Then we adopt a Natural
Gradient Boost (NGBoost) method to analyze the relationship between the SHARP
parameters and the successive flare bursts. Based on the correlation analysis
and the importance distribution returned from NGBoost, we select 8 most flare-
related SHARP parameters. Finally, we discuss the physical meanings of the 8
selected parameters and their relationship with flare production.

##  1 Introduction

As one of the most prominent eruptive phenomena, solar flares can produce
significant space weather effects including energetic particles. In flare-
eruptive active regions (ARs), flares may occur successively. Since solar
flares are often accompanied by coronal mass ejections (CMEs), successive
flare bursts may appear with successive CMEs, the interaction between which
may result in extreme space weather (Liu et al., 2014b, 2019). Establishing a
practical flare forecast system, especially for successive flare bursts in the
same AR, is thus of importance for space weather. Since the vector magnetic
field can only be mapped directly in the photosphere, using photospheric
magnetic field data to parameterize ARs has been widely adopted in studies
about flare prediction(e.g., Yang et al., 2013; Bobra & Couvidat, 2015; Wang
et al., 2019, 2020; Cicogna et al., 2021; Yi et al., 2021; Ribeiro & Gradvohl,
2021). Due to the complex current distribution (Török et al., 2014) and the
magnetic field evolution, it is difficult to derive a single parameter that
governs flare production in the same AR. Therefore, identifying parameters
that are the most relevant to flare productivity is of great importance.

The _Helioseismic and Magnetic Imager_ (HMI; Schou et al., 2012), aboard the
_Solar Dynamics Observatory_ (SDO), continuously maps the full-disk
photospheric vector magnetic field since 2010 May 1 every 12 minutes. It also
provides a data product named Space-Weather HMI Active Region Patch (SHARP),
which offers a series of key parameters derived from the magnetic field in ARs
(Bobra et al., 2014; Hoeksema et al., 2014). Lee et al. (2018) applied
$I_{8d}$, the sum of the GOES soft X-ray peak fluxes of flares in the same AR
for 8 days, to investigate the correlation between SHARP parameters and flare
productivity in time-series, concluding that the SHARP parameters are to some
extent related to flare productivity. Wang et al. (2019) calculated the SHARP
parameters in the polarity inversion line area of an AR, and employed a
Support Vector Machine (SVM) algorithm to perform flare prediction with the
modified SHARP parameters. Many of the previous flare prediction studies focus
on determining if an AR is flare-eruptive or flare-quiet, instead of looking
into the successive flare bursts along with the evolution of the AR (e.g.,
Yang et al., 2013; Bobra & Couvidat, 2015; Wang et al., 2019; Cicogna et al.,
2021; Yi et al., 2021; Ribeiro & Gradvohl, 2021). Some studies considered the
successive flare production, but focused on only limited parameters such as
the degree of current neutralization (e.g., Kontogiannis et al., 2017). Here
our approach is different. We consider all the flares above C1.0 in the same
AR, and study the relationship between the successive flare bursts and SHARP
parameters.

The polarity inversion line (PIL), where the polarity of the magnetic field is
reversed, has been proved to have strong relationships with solar flares and
CMEs (e.g., Falconer et al., 2002; Mason & Hoeksema, 2010; Welsch et al.,
2011; Georgoulis et al., 2012; Moore et al., 2012; Török et al., 2014; Liu et
al., 2017a; Wang et al., 2018, 2022). Parameters, such as the shear angle and
the emerging magnetic flux, have been proved to have a proportional
relationship with flare productivity, if narrowed to the PIL area (Moore et
al., 1987; Schrijver, 2007). Previous works suggest that the PIL modified
SHARP parameters perform better than original SHARP parameters with respect to
flare prediction(e.g., Wang et al., 2019; Cicogna et al., 2021). In our study,
we also employ similar methods to see the effect of PIL modification when
studying the successive flares in the same ARs.

In this paper, we evaluate the relationship between SHARP parameters and
successive flare production in two ways. The first method is the calculation
of correlation coefficients. The second method is based on a machine learning
algorithm named Natural Gradient Boost (NGBoost), which is a gradient boost
algorithm that can be easily used without time-consuming hyperparameter tuning
(Duan et al., 2020). This algorithm returns a probability distribution for the
target and an importance distribution of the input parameters. In this study,
the input parameters are the SHARP parameters. With the correlation result and
the importance distribution returned from NGBoost, we could identify the most
flare-related SHARP parameters, and make a physical analysis.

The paper is organized as follows. In Section 2, we describe the data and the
methodology that are employed in this paper. In Section 3, we show the
correlation results, the importance distribution of the SHARP parameters, and
the test of the selection. The underlying physics analysis of the selected
parameters is also discussed in Section 3. The conclusion is given in Section
4.

##  2 Data and Methodology

###  2.1 Active Region Selection

We have selected 6 ARs, which all have at least one major flare above level
X1.0. Table 1 lists the index, start and end time, number of flares, maximum
flare, and the number of groups of SHARP data that we have employed for the
ARs.

Table 1: List of active regions in our study. AR | Start time | End time | Number of flares | Maximum flare | Groups of data
---|---|---|---|---|---
11158 | 2011-02-13 11:58:12 | 2011-02-18 11:46:12 | 41 | X2.2 | 600
11166 | 2011-03-08 23:58:14 | 2011-03-11 23:46:15 | 12 | X1.5 | 357
11283 | 2011-09-05 23:58:22 | 2011-09-10 23:46:21 | 13 | X2.1 | 599
11429 | 2012-03-04 15:58:14 | 2012-03-12 13:58:15 | 45 | X5.4 | 900
12017 | 2014-03-26 23:58:16 | 2014-03-31 23:46:17 | 15 | X1.0 | 588
12673 | 2017-09-01 23:58:42 | 2017-09-07 18:46:41 | 42 | X9.3 | 620

Since the ARs are at different latitudes, and the longitudes of them are also
different when they are formed, the time durations for them to stay on the
visible side of the Sun are different, resulting in differences in the
effective observation times.

The ARs have been studied widely from many aspects, such as the energy release
process (e.g., Tarr et al., 2013; Feng et al., 2013; Wang et al., 2014, 2018,
2022), the structure of the nonlinear force-free field (e.g., Feng et al.,
2013; Wang et al., 2013; Vemareddy & Wiegelmann, 2014), and the CMEs
associated with the major flares (e.g., Liu et al., 2013, 2014a).

###  2.2 Flare “Envelope” Curve

The intensity of flares is defined by GOES Soft X-ray (SXR) flux, and we only
select the flares bigger than C1.0 that has a peak flux of $1\times 10^{-6}\
W/m^{2}$. We ignore the temporal evolution of an individual flare, and define
its time tag as the flare peak time. The curve “enveloping” successive flares
is calculated from:

$f_{i}=log_{10}(\Phi_{i})+6,$ (1)

where $f_{i}$ is the defined flare index in the flare “envelope”, and
$\Phi_{i}$ is the corresponding peak flux. The flare data and the GOES SXR
data are obtained from the SolarflareDatabase
111https://solarflare.njit.edu/datasources.html. An example of the flare
“envelope” of AR 12673 is shown in Figure 1b.

Figure 1: The original GOES SXR curve and the flare “envelope” of AR 12673.
Panel “a” refers to the original GOES SXR curve, and panel “b” represents the
corresponding flare “envelope”. The temporal resolution of both panels is
720s, which is the same as that of the SHARP data.

The use of the flare envelope allows us to derive the flare data in a specific
AR instead of the whole solar disk, and set a threshold for the flares
conveniently. Flares from other ARs can be readily eliminated. It also allows
us to easily adjust the temporal resolution of the flare curve to match the
temporal resolution of SHARP. The main disadvantage of the envelope curve is
that it is obtained from discrete flare bursts. Temporal details of individual
flares are ignored when the flare “envelope” curve is employed.

Note that the flare envelope curve is less zigzag than the curves derived from
the original GOES SXR flux curve. For example, we have tried another curve
that is directly derived from the GOES SXR flux curve (see Figure 1a). The
GOES SXR data points are kept if they are within the timespan of the flares
from the selected AR. The removed data points are filled by a linear
interpolation. Compared to the flare envelope curve, the new curve is not
significantly different in shape but much more zigzag. A more zigzag curve of
flares will increase the difficulty for correlation analysis and further
NGBoost analysis.

###  2.3 Space-Weather HMI Active Region Patch (SHARP)

SHARP222http://jsoc.standford.edu/ajax/lookdata.html provides a series of key
parameters derived from photospheric magnetic fields. We have selected 16
typical SHARP parameters (Bobra et al., 2014). Table 2 lists the keywords, the
descriptions and the calculation formulas for all the SHARP parameters we have
used in this study. The parameters are generally adopted from Leka & Barnes
(2003a).

Table 2: Brief descriptions and corresponding formulas of the SHARP parameters (Bobra et al., 2014). Keyword | Description | Calculation formula
---|---|---
TOTUSJH | Total unsigned current helicity | $\bm{H}_{c_{total}}\propto\sum|\bm{B_{z}}\cdot\bm{J_{z}}|$
TOTPOT | Total photospheric magnetic free energy density | $\rho_{tot}\propto\sum(\bm{B}^{obs}-\bm{B}^{pot})^{2}dA$
TOTUSJZ | Total unsigned vertical current | $\bm{J}_{z_{total}}=\sum|\bm{J}_{z}|dA$
ABSNJZH | Absolute value of net current helicity | $\bm{H}_{c_{abs}}\propto|\sum\bm{B}_{z}\cdot\bm{J}_{z}|$
SAVNCPP | Sum of modulus of the net current per polarity | $\bm{J}_{z_{sum}}\propto|\sum_{B_{z}^{+}}\bm{J}_{z}dA|+|\sum_{B_{z}^{-}}\bm{J}_{z}dA|$
USFLUX | Total unsigned flux | $\Phi=\sum|\bm{B}_{z}dA|$
MEANPOT | Mean photospheric magnetic free energy | $\overline{\rho}\propto\frac{1}{N}\sum(\bm{B}^{obs}-\bm{B}^{pot})^{2}$
R_VALUE | Sum of flux near PIL |  $\Phi=\sum|\bm{B}_{z}|dA$ within R mask
MEANSHR | Mean shear angle | $\overline{\Gamma}=\frac{1}{N}\sum cos^{-1}(\frac{\bm{B}^{obs}\cdot\bm{B}^{pot}}{|\bm{B}^{obs}||\bm{B}^{pot}|})$
MEANGAM | Mean angle of field from radial | $\overline{\gamma}=\frac{1}{N}\sum tan^{-1}\frac{B_{h}}{B_{z}}$
MEANGBT | Mean gradient of total field | $\overline{|\nabla\bm{B}_{tot}|}=\frac{1}{N}\sum[(\frac{\partial B_{tot}}{\partial x})^{2}+(\frac{\partial B_{tot}}{\partial y})^{2}]^{(1/2)}$
MEANGBZ | Mean gradient of vertical field | $\overline{|\nabla\bm{B}_{z}|}=\frac{1}{N}\sum[(\frac{\partial B_{z}}{\partial x})^{2}+(\frac{\partial B_{z}}{\partial y})^{2}]^{(1/2)}$
MEANGBH | Mean gradient of horizontal field | $\overline{|\nabla\bm{B}_{h}|}=\frac{1}{N}\sum[(\frac{\partial B_{h}}{\partial x})^{2}+(\frac{\partial B_{h}}{\partial y})^{2}]^{(1/2)}$
MEANJZH | Mean current helicity ($B_{z}$ contribution) | $\overline{\bm{H}_{c}}\propto\frac{1}{N}\sum\bm{B}_{z}\cdot\bm{J}_{z}$
MEANJZD | Mean vertical current density | $\overline{\bm{J}_{z}}\propto\frac{1}{N}\sum(\frac{\partial B_{y}}{\partial x}-\frac{\partial B_{x}}{\partial y})$
MEANALP | Mean characteristic twist parameter: $\alpha$ | $\overline{\alpha}\propto\frac{\sum\bm{J}_{z}\cdot\bm{B}_{z}}{\sum\bm{B}_{z}^{2}}$

###  2.4 Correlation and NGBoost Analysis

First, we look at the SHARP parameters for the whole AR, and calculate the
correlation coefficients between the SHARP curves and the flare envelope. We
also look at the SHARP parameters in the PIL area to study how PIL
modification affects the correlation results. The method to identify the PIL
initially comes from Schrijver (2007). To be specific, we set a threshold of
$2\sigma$ (200G) (Hoeksema et al., 2014), any positive field pixel with
$B_{r}>200$G is set to 1 and the rest are 0 to make the positive bitmap, and
any negative field pixel with $B_{r}<-200$G is set to -1 and the rest are 0 to
make the negative bitmap. Then we adopt a cluster algorithm called Density-
Based Spatial Clustering of Applications with Noise (DBSCAN; Ester et al.,
1996) to remove the noise data points and leave only the main parts of the
positive and negative bitmaps. We finally multiply the two bitmaps after
convolving them with a Gaussian kernel to generate the PIL mask.

We also adopt a Natural Gradient Boost (NGBoost)
method333https://stanfordmlgroup.github.io/projects/ngboost (Duan et al.,
2020) to analyze the relationship between flare production and SHARP
parameters and obtain the importance distribution. The NGBoost method is a
gradient boosting algorithm for probabilistic regression problems, which
returns a possibility distribution of the target. The possibility distribution
is defined by two parameters, “location” and “scale”. The “location” parameter
represents the location of the maximum point of the distribution. It shifts
the graph, relative to the normal probability distribution function (PDF),
left or right on the horizontal axis depending on if its value is negative or
positive. The “scale” parameter defines the width of the distribution. It
expands or shrinks the graph, relative to the normal PDF, depending on whether
its value is larger or smaller than 1.

We apply the NGBoost method to SHARP parameters as follows. We discard the
first 2 hours of flare envelope data, and match the remaining flare envelope
data to the original SHARP parameters, in order to establish a 2-hour time
difference between the flare envelope and the SHARP curve. This is to explore
whether the SHARP parameters before the flares can reflect flare productivity.
Then, we take the first 10 groups of data, which includes the first 10 groups
of SHARP parameters and the first 10 flare envelope points, to be the initial
training set. We start the calculation from the 11th groups of the SHARP
parameters. With the training of the first 10 groups of data finished, we
input the feature vector of the 11th group of SHARP parameters to the trained
model, and obtain the corresponding “flare envelope” point calculated from
NGBoost. An importance distribution of the SHARP parameters in this training
step is returned from NGBoost as well. In the next step, we take the first 11
groups of data as the training set, input the 12th SHARP parameters to the
trained model, and obtain the next calculated “flare envelope” point along
with a new importance distribution. This operation is conducted to the last
group of data. Finally, we attain a new “flare envelope” calculated from
NGBoost and the importance distributions for all the steps. The final
importance distribution of each SHARP parameter is the average value of all
the importance distributions that we have obtained. The explanation of the
NGBoost method can be found in the Appendix. It is worth mentioning that the
2-hour difference between the flare envelope and the SHARP parameters is only
used in the “NGBoost Analysis” section (Section 3.2). In the other two
sections of Section 3, the SHARP parameters and the flare envelope are
simultaneous.

##  3 Results

###  3.1 Correlation Analysis

As mentioned above, we first look at the SHARP parameters in the whole AR and
calculate their correlation coefficients with the flare envelope to estimate
how each SHARP parameter is related to the flare envelope. We take the AR
12673 as an example in Figure 2. As can be seen in Figure 2, in AR 12673, the
range of the correlation coefficients (CCs) is wide for different SHARP
parameters. For example, the CC for “TOTPOT” is over 0.8, while the absolute
value of the CC for “MEANGBH” is only 0.104, which suggests that “MEANGBH” is
much less relevant to the flare envelope than “TOTPOT” in this AR.

Figure 2: The 16 SHARP parameters and the flare envelope for AR 12673. The
black curve shows different SHARP parameters and the green curve is the flare
envelope. The left y-axises represent the values of the SHARP parameters, and
the right y-axises represent the flare index. The correlation coefficients
(CCs) are shown at the top of each panel. The SHARP parameters are normalized.

We calculate the CCs for all the SHARP parameters, in the whole AR and in the
PIL area, for all the ARs that we have selected. Then we plot the average and
variance of the CCs for SHARP parameters in Figure 3. Figure 3 consists of two
panels, referring to the SHARP parameters in the whole AR and the SHARP
parameters in the PIL area. We mark the parameters that we think have the best
performance in correlation with the blue ellipse in the left panel, which
includes “MEANJZH”, “MEANGBH”, “MEANALP”, “USFLUX”, “MEANPOT”, “MEANSHR” and
“MEANGAM”.

Figure 3: Scatter plot of the average and variance of the correlation
coefficients between the SHARP parameters and the flare envelope. The left
panel refers to the whole AR case, and the right panel refers to the PIL-
modified case. The x-axis represents the average value, and y-axis represents
the variance value. A parameter has a better performance if it is closer to
the lower right corner of the figure. Inside the blue ellipse are the
parameters that we think have the best performance in correlation.

Our results in Section 3.3 present that, in agreement with previous studies
(e.g., Schrijver, 2007; Wang et al., 2019), the SHARP parameters that are
believed to be more related to successive flare occurrence do show more
significant changes in the PIL area when major flares happen. However, the two
panels in Figure 3 suggest that the average values of some SHARP parameters in
the PIL area, such as “MEANPOT” and “MEANJZH”, are lower than those in the
whole AR. This is due to sharper fluctuations of the SHARP curves in the PIL
area. The value of the correlation coefficient between the two curves is not
only affected by their overall shapes, but also their smoothness. Therefore,
the SHARP curves in the PIL area may fit the flare envelope better in the
overall shape, but the correlation coefficients may be lower due to their
stronger oscillation. The more intense fluctuations of the SHARP curves in the
PIL area can be attributed to two reasons: (i) During the evolution of a big
AR, the coverage of the PIL area varies, which results in fluctuations in
SHARP parameters calculated in the PIL area. (ii) Some ARs have more than one
PIL at one time, while some other ARs may have no obvious PIL. Even in the
same AR, such as AR 12673, the quantity of PILs may change with evolution.
Figure 4 presents the evolution of the PILs in AR 12673. AR 12673 starts with
no obvious PIL (see Figure 4a). However, as it becomes more active, it has
more than one PIL and the coverage of the PILs increases significantly (see
Figure 4 (b-c)). Panel “d” in Figure 4 also shows that the coverage of the
PILs fluctuates strongly when the AR gets complicated, which may cause severe
fluctuations in the SHARP curve. Similar features about the change in the
coverage and the number of PILs can be found in all the ARs that we select.

Figure 4: Illustration of the change in the coverage and the number of PILs in
AR 12673. Panel (a), panel (b) and panel (c) are the composite maps of the Br
map and their PILs at different times. The PIL area is marked by the purple
lines. Panel (d) shows the change in the number of pixels of PILs in AR 12673.
The three vertical lines in panel (d) correspond to the times of panels (a-c).

###  3.2 NGBoost Analysis

Since the correlation coefficients between the flare envelope and the SHARP in
the whole AR are usually higher, we only apply the NGBoost method to the SHARP
parameters in the whole AR. The result is shown in Figure 5. The blue curve is
the “flare envelope” we calculated from NGBoost, and the green curve is the
original flare envelope. The NGBoost curves show good fits to the flare
envelopes in all the ARs we select.

Figure 5: The original flare envelope and the corresponding “flare envelope”
calculated from NGBoost for the 6 selected ARs. The green curves are the
original flare envelopes, and the blue curves are from NGBoost.

With the importance distributions of the SHARP parameters returned from
NGBoost, we calculate the ranking of different SHARP parameters as shown in
Figure 6. It contains the ranking average and corresponding variance of all
the SHARP parameters in the 6 ARs. The x-axis stands for the ranking average
and the y-axis stands for the ranking variance, both of which are better when
smaller. Thus, if a parameter is closer to the bottom left corner, it is of
more importance. As Figure 6 suggests, we can derive “TOTPOT”, “MEANPOT”,
“USFLUX” and “MEANGAM” out of the 16 parameters as the good performance group.
Combining the correlation analysis in Figure 3, we put another 4 SHARP
parameters into the good performance group. The 4 SHARP parameters are
“MEANJZH”, “MEANGBH”, “MEANALP”, “MEANSHR”. All the selected parameters have
been marked by the blue ellipses.

Before we start to analyze the 8 parameters separately, we have verified the
reliability of the 8 parameters by applying the NGBoost method once more to
the 8 parameters instead of all the SHARP parameters. The result is shown in
Figure 7. We have also conducted the same operation with fewer SHARP
parameters, and the fitting of the two curves is not as good as the eight-
element result. This suggests that the selection is successful and convincing.

###  3.3 Underlying physics of the selected SHARP parameters

Figure 6: Scatter diagrams for the ranking of different SHARP parameters
calculated from NGBoost. The left panel shows the ranking for the parameter
“location” in the possibility distribution, and the right one shows the
ranking for the parameter “scale”. Inside the blue ellipse are the parameters
of the greatest importance. Figure 7: Comparison between the original flare
envelope and the “flare envelope” curve calculated from NGBoost with the
selected 8 SHARP parameters only.

We analyze the 8 selected SHARP parameters in this section. Figure 8 shows an
example for the evolution of the selected SHARP parameters in the whole AR and
the PIL area for AR 11283. Two adjacent panels, such as the panel “a” and
panel “b”, represent the same parameter. The left panel shows the SHARP in the
whole AR, and the right one shows the SHARP in the PIL area. Here we include
the PIL-modified SHARP parameters to study how the SHARP parameters evolve in
the PIL area during flares.

The first two parameters are “TOTPOT” and “MEANPOT”, which have similar
physical meanings but are derived in different ways. “TOTPOT” is the total
value of the magnetic free energy, and “MEANPOT” is the average value. Free
energy is highly related to the nonlinear-force-free field. According to
previous studies, a great amount of accumulated magnetic free energy is
released within a short time to provide the energy for flares and CMEs (e.g.,
Antiochos et al., 1999; Amari et al., 2000; Chen & Shibata, 2000; Wang et al.,
2014). Guo et al. (2008) analyzed the flare in December 2006 by extrapolating
its 3D magnetic field configuration. They calculated the free energy and the
shear angle of the AR, studied the energy release process throughout the
evolution of the magnetic field, and explained the energy release process
during a flare. Schrijver et al. (2008) also explained the free energy release
process by extrapolating the nonlinear-force-free field, and found that a free
energy drop of $\sim 10^{32}$ erg happened as the X flare occurred. In our
study, the two SHARP parameters about the free energy are also considered as
two of the most related parameters for flare occurrence.

In Figure 8, panel “a” and panel “b” represent “TOTPOT”, and panel “c” and
panel “d” represent “MEANPOT”. We can see in the panels that, comparing to the
two parameters in the whole AR, “TOTPOT” and “MEANPOT” show a more obvious
increase in the PIL area when big flares happen. The absolute value of
“MEANPOT” is much higher than that in the whole AR. These features are found
in all the ARs we select, which suggests that the storing and releasing
process of free energy is more obvious in the PIL area. This conclusion is in
agreement with previous studies (e.g., Cicogna et al., 2021; Wang et al.,
2019).

The next parameter that we look into is “USFLUX”, which represents the total
unsigned flux. “USFLUX” is shown in panels (e-f) in Figure 8. The correlation
between the total unsigned flux and flare productivity has been revealed by
previous studies (e.g., Schrijver, 2007; Vasantharaju et al., 2018; Bobra &
Couvidat, 2015; Kazachenko et al., 2017). Vasantharaju et al. (2018) made a
statistical study of some magnetic imprints and their correspondence with
flare productivity. Following the results of previous studies (e.g.,
Schrijver, 2007; Kazachenko et al., 2017), they took three kinds of total
unsigned magnetic flux into consideration, the flux of the whole AR, the flux
along the PIL, and the flux in the flare ribbon area. Comparing to the total
unsigned flux in the whole AR, a slight and a strong increase in correlation
with flares are found along with the PIL and in the flare ribbon area,
respectively.

Figure 8: The 8 parameters in the whole AR and in the PIL area for AR 11283.
Two adjacent panels, such as panel “a” and “b”, belong to a group of the same
SHARP. In the same group, the left panel is the SHARP in the whole AR, and the
right one is the SHARP in the PIL area. Counting from left to right and up to
down, the 8 groups are: “TOTPOT”, “MEANPOT”, “USFLUX”, “MEANGAM”, “MEANGBH”,
“MEANJZH”, “MEANSHR” and “MEANALP”. Each panel is labeled with a letter from
“a” to “p”. The SHARP parameters are normalized. The corresponding correlation
coefficient (CC) is at the top of each panel. The maximum value of each SHARP
parameter is also at the top. Similar format to Figure 2.

The two parameters of “MEANGAM” and “MEANGBH” are both related to the
horizontal magnetic field. In Figure 8, panel “g” and panel “h” represent
“MEANGAM”, and panel “i” and panel “j” represent “MEANGBH”. The main finding
about the horizontal magnetic field is that the horizontal field tends to
increase significantly in the PIL area after flare eruptions, indicating that
the vector magnetic field tends to be more horizontal after flares (Petrie,
2012a, b; Sun et al., 2012; Wang et al., 2014). This is consistent with our
conclusion drawn from Figure 8. Compared to panel “g” and panel “i”, the SHARP
curves in panel “h” and panel “j” show more obvious increases in the PIL area
when big flares happen. This feature of these two parameters is also found in
all the ARs we select. Previous studies suggest that a rapid and irreversible
increase appears in the horizontal magnetic field after an eruption, and it
becomes more aligned with the PIL as well. Song & Zhang (2016) performed a
detailed study by considering the continuum intensity. They compared the
continuum intensity and the horizontal magnetic field in 5 ARs, and revealed
that the differences of the continuum intensity and the horizontal magnetic
field before and after the major flares (above X1.0) are strongly correlated
both spatially and temporally in the outer penumbra of the sunspot. A similar
conclusion has been drawn by Lu et al. (2019) through a statistical analysis
of 20 X-class flares. All these works indicate that the horizontal magnetic
field is of great importance for successive flares in the same AR.

The parameter ”MEANJZH”, shown in panels (k-l) in Figure 8, stands for the
mean current helicity with $B_{z}$ contribution. Current helicity is crucial
to magnetic energy build-up. An explanation for the relationship between
flares and current helicity is the “$\alpha-$effect”, indicating that the
energy of small-scale fluctuations of magnetic and velocity fields could be
transferred into the energy of large-scale currents (Seehafer, 1990). The
necessary condition for such effect may be the presence of the predominant
sign of the electric helicity over the whole AR. Abramenko et al. (1996)
conducted a statistical study involving 40 ARs and suggested that the
condition of the “$\alpha-$effect” worked in over $90\%$ of all the cases.
Yeates et al. (2008) performed a simulation to study the current helicity.
They discussed the sources of the current helicity and suggested that field
lines are skewed as they cross the PIL. Similar to current helicity, magnetic
helicity is also believed to be closely related to the mechanisms that drive
solar activities. “MEANALP”, shown in panels “o” to “p” in Figure 8,
represents the mean characteristic twist parameter. The method used to
calculate this SHARP parameter was introduced by Hagino & Sakurai (2004) and
adopted by Bobra et al. (2014) in the SHARP package. Wang et al. (2022)
studied the 4 homologous events (all are eruptive flares) in AR 11283. They
analyzed the buildup process of the magnetic flux ropes (MFRs), especially the
evolution of the magnetic helicity and the energy during the 4 eruptions.
Their results show that the MFRs, which are built up along the collisional
polarity inversion lines (cPILs), store a large amount of energy and magnetic
helicity. The buildup process is involved with magnetic cancellation along the
cPIL via collisional shearing, which may contribute to solar eruptions
(Chintzoglou et al., 2019). Li et al. (2022) proposed several new parameters
and investigated their ability in distinguishing eruptive and confined solar
flares. They found that the ratio between the “MEANALP” in the flaring PILs
and the total unsigned flux in the whole AR ($\alpha_{FPIL}/\Phi_{AR}$) is the
most capable parameter to distinguish ARs that can produce eruptive flares. In
our work, as shown in Figure 8, the values of the parameters “MEANJZH” and
“MEANALP” in the PIL area are about ten times the values in the whole AR,
indicating that both current and magnetic helicities accumulate in the PIL
area. This is in agreement with the mentioned studies.

The last parameter we select is “MEANSHR”, which stands for the mean shear
angle. The shear angle is defined as the angle between the potential field and
the observed field in the photosphere. It is commonly accepted that a more
strongly sheared magnetic filed stores more free energy. Kusano et al. (2012)
conducted a simulation work to interpret the types of magnetic field
structures that can drive solar eruptions. Their simulation results show that
the shear angles of large-scale eruptions accurately indicate the kinetic
energy produced by eruptions. More sheared magnetic fields favor the onset of
large-scale eruptions. A lot of studies have been published to reveal the
relationship between shear angles and flare production, suggesting that the
shear angle is important in studying flare bursts (e.g., Moore et al., 1987;
Zhang, 2001; Leka & Barnes, 2007; Guo et al., 2008; Wang et al., 2022; Li et
al., 2022).

##  4 Summary and Conclusions

This paper has identified the SHARP parameters that are most related to
successive flare bursts in the same AR. For this purpose, we adopt the
envelope of the successive flare bursts to represent discrete individual
events, and calculate the correlation coefficients between the envelope and
the SHARP parameters. The NGBoost method (Duan et al., 2020) is also employed
to obtain the importance distribution of the SHARP parameters. With the
importance distribution from NGBoost and the correlation results, we have
selected 8 parameters. A time difference of 2 hours between the SHARP
parameters and the flare envelope is established in the NGBoost section, which
is for exploring if the SHARP parameters before flares could reflect flare
productivity. In fact, we have conducted the NGBoost analysis 10 times with
time differences varying from 12 minutes to 2 hours, and the selections of
SHARP parameters are almost the same. By comparing our results with previous
works, we analyze the physical meanings of the 8 parameters and their
relationships with flare productivity. The major conclusions are summarized
below.

  1. 1.

The most relevant SHARP parameters to the successive flare productivity in the
same AR are “TOTPOT”, “MEANPOT”, “USFLUX”, “MEANGAM”, “MEANGBH”, “MEANJZH”,
“MEANSHR”, and “MENAALP”. These parameters are related to the free energy,
unsigned magnetic flux, horizontal magnetic fields, current and magnetic
helicity, and shear angles, all of which are confirmed to have strong
connections with flares.

  2. 2.

The PIL is important in studying successive flares in the same AR. We do see
more significant changes in the selected parameters in the PIL area, in
particular during big flares. A significant store and release process of the
free energy before and after flares can be found in both the PIL area and the
whole AR, but it is more concentrated in the PIL area. The horizontal field
becomes more dominant after flare eruptions in the PIL area, and the magnetic
field is also more twisted in the PIL area. However, during the evolution of
an AR, the coverage and number of PILs may vary significantly. This would
increase the difficulty of investigating the mechanisms that drive successive
solar eruptions.

  3. 3.

The establishment of a practical flare forecast system requires full
understanding of the physical mechanism of flare bursts. A time difference
between the parameters and the flare indexes may help us to explore the
trigger of flares. We have made an example with time difference of two hours,
which seems promising, but more work is needed in the future.

A lot of work has been put into identifying the parameters that carry the most
information for flare bursts in order to lay foundation for a practical
prediction scheme. Leka & Barnes (2007) conducted a study on distinguishing
the photospheric magnetic field properties that are important for flare
production. Their work shows that the properties related to free energy, total
unsigned flux, vertical currents and shear angles are strongly related to
flare production. Bobra & Couvidat (2015) found that the magnetic energy,
vertical current and current helicity are the most useful physical quantities
in flare prediction, in agreement with Leka & Barnes (2003a) and our result
here. Ahmed et al. (2013) employed a Solar Monitor Active Region Tracker
(SMART) to estimate the efficiency of such parameters in flare prediction, and
selected 6 physical quantities, which are believed to carry the most
information for flare bursting, out of 21 parameters. Unlike the previous
studies, we focus on successive flare bursts in the same AR, and identify the
most flare-related SHARP parameters with a time difference of 2 hours. It
still remains a problem to give a practical flare forecasts in the
continuously evolving ARs with a high accuracy. Interactions between adjacent
flares have highly increased the difficulty of establishing a precise
prediction method. More work about the physical mechanisms of the flares is
needed to establish a truly reliable flare forecast system.

## Acknowledgements

This research was supported by the National Key R$\&$D Program of China
No.2021YFA0718600, NSFC (41774179, 12073032, 42004145, 42150105), the
Strategic Priority Research Program of Chinese Academy of Sciences
(XDA15018500), and the Specialized Research Fund for State Key Laboratories.
Y.G. was supported by NSFC (11773016, 11961131002, and 11533005) and
2020YFC2201201. H.R. acknowledges helpful discussions with Dr. Huidong Hu. We
acknowledge the work of the Sunpy team (The SunPy Community et al., 2020). We
are grateful to the work of Glogowski et al. (2019), whose Python package
_drms_ has offered us a convenient tool to access the SDO data. The data used
here are courtesy of NASA/SDO and the HMI science team, as well as the
Geostationary Satellite System (_GOES_) team.

## The NGBoost method

NGBoost is a gradient boosting algorithm that aims at estimating the
parameters of a probability distribution $P_{\theta}(y|\bm{x})$, where
$\bm{x}$ is the feature vector of observation, and $y$ is the target variable.
In our work, $\bm{x}$ is the vector of SHARP parameters and $y$ is the
corresponding flare envelope value. NGBoost can also be applied to point
estimation by extracting point expectations $\widehat{\mathbb{E}}[y|x]$ from
the estimated distribution $\widehat{P}_{\theta}(y|x)$ (Duan et al., 2020).
There are three main modules in NGBoost: (i) the base learner $f$, (ii) the
probability distribution $P_{\theta}$, and (iii) a proper scoring rule. Figure
2 in Duan et al. (2020) shows the schematic illustration of NGBoost. The base
learners are crucial in gradient boosting algorithms. In general, gradient
boosting algorithms rely on the sequential training of the base learners to
form an additive ensemble. The ensemble of the previous learners estimates a
current residual, which will then be applied to optimize the new learner. The
output of that learner is scaled by a learning rate, and it is then appended
to the current ensemble. In our work, we employ the decision tree that is
suggested in Duan et al. (2020) as the base learner.

Regarding the probability distribution, we choose the normal probability
distribution. It is defined by $\theta$, which is a vector composed of
“location” and “scale”. The next step is to determine a proper scoring rule. A
proper scoring rule generates a rating between the estimated distribution and
the observation, then returns the rating to the base learner such that the
distribution of the outcomes gets the best score in expectation. The most
commonly used scoring rule, which is also applied in this context, is the
logarithmic score:

$\mathcal{L}(\theta,y)=-log(P_{\theta}(y)).$ (1)

Another important part of NGBoost is the natural gradient
$\tilde{\bigtriangledown}_{\theta}$. The natural gradient is motivated by
information geometry. It makes the optimization problem invariant to the
parameterization, and it also makes the learning more stable and effective
(Amari, 1998).

Figure 9: The importance feature of the SHARP parameters in AR 12673. The left
panel refers to the importance feature distribution of “location”, and the
right panel refers to the importance feature distribution of “scale”.

To be specific, assuming that we have a dataset
$D=\\{x_{i},y_{i}\\}_{i=1}^{M}$, where $M$ is the number of the training
examples. The algorithm first establishes an initial predicted parameter
$\theta^{(0)}$, which is defined as:

$\theta^{(0)}=argmin_{\theta}\sum_{i=1}^{M}\mathcal{L}(\theta,y_{i}).$ (2)

The initial predicted parameter is a standard set of parameters that aims at
minimizing the sum of the scoring rule $\mathcal{L}$. Next, for each iteration
n, we estimate the individual natural gradient for each example i as:

$g_{i}^{(n)}=\mathcal{I}_{\mathcal{L}}(\theta_{i}^{(n-1)})^{-1}\cdot\bigtriangledown_{\theta}\mathcal{L}(\theta_{i}^{(n-1)},y_{i}),$
(3)

where the $\mathcal{I}_{\mathcal{L}}(\theta)$ is the Fisher Information
carried by an observation about $P_{\theta}$. For logarithmic score, the
Fisher Information is defined as:

$\mathcal{I}_{\mathcal{L}}(\theta)=\mathbb{E}_{y\thicksim
P_{\theta}}[\bigtriangledown_{\theta}\mathcal{L}(\theta,y)\bigtriangledown_{\theta}\mathcal{L}(\theta,y)^{T}].$
(4)

After being established, the natural gradient and the input vectors
$\boldsymbol{x}_{i}$ are used to train the base learners of iteration
$f^{(n)}$. The training process returns a fitted base learner that is the
projection of the natural gradient. A scale factor $\rho^{(n)}$ is needed to
scale the projected gradient. It is computed as:

$\rho^{(n)}=argmin_{\rho}\sum_{i=1}^{M}\mathcal{L}(\theta_{i}^{(n-1)}-\rho\cdot
f^{(n)}(\bm{x}_{i}),y_{i}),$ (5)

where $\rho$ starts from $\rho=1$. Then the parameters $\theta_{i}^{(n)}$ are
updated by:

$\theta_{i}^{(n)}=\theta_{i}^{(n-1)}-\eta(\rho^{(n)}f^{(n)}(\bm{x}_{i})),$ (6)

where $\eta$ denotes a usual learning rate (typically 0.1 or 0.01).

The use of NGBoost in this context is explained in Section 2.4. We focus on
the importance feature of each SHARP parameter. The SHARP parameter that has
higher importance feature value in the NGBoost analysis is believed to be more
responsible for successive flare occurrence. An example of the importance
feature distribution is shown in Figure 9. The NGBoost method can be
conveniently used and interpreted
444https://stanfordmlgroup.github.io/ngboost/intro.html with the python
package _ngboost_ 555https://github.com/stanfordmlgroup/ngboost.

## References

  * Abramenko et al. (1996) Abramenko, V. I., Wang, T., & Yurchishin, V. B. 1996, SoPh, 168, 75
  * Ahmed et al. (2013) Ahmed, O. W., Qahwaji, R., Colak, T., et al. 2013, SoPh, 283, 157
  * Amari (1998) Amari, S.-I. 1998, Neural computation, 10, 251
  * Amari et al. (2000) Amari, T., Luciani, J. F., Mikic, Z., & Linker, J. 2000, ApJ, 529, L49
  * Antiochos et al. (1999) Antiochos, S. K., DeVore, C. R., & Klimchuk, J. A. 1999, ApJ, 510, 485
  * Bobra & Couvidat (2015) Bobra, M. G., & Couvidat, S. 2015, ApJ, 798, 135, doi: 10.1088/0004-637x/798/2/135
  * Bobra et al. (2014) Bobra, M. G., Sun, X., Hoeksema, J. T., et al. 2014, SoPh, 289, 3549, doi: 10.1007/s11207-014-0529-3
  * Chen & Shibata (2000) Chen, P. F., & Shibata, K. 2000, ApJ, 545, 524
  * Chintzoglou et al. (2019) Chintzoglou, G., Zhang, J., Cheung, M. C. M., & Kazachenko, M. 2019, ApJ, 871, 67
  * Cicogna et al. (2021) Cicogna, D., Berrilli, F., Calchetti, D., et al. 2021, arXiv preprint arXiv:2105.00897
  * Duan et al. (2020) Duan, T., Avati, A., Ding, D. Y., et al. 2020, PMLR, 2690
  * Ester et al. (1996) Ester, M., Kriegel, H.-P., Sander, J., & Xu, X. 1996, Int. Conf. Knowledge Discovery and Data Mining, 240, 6
  * Falconer et al. (2002) Falconer, D. A., Moore, R. L., & Gary, G. A. 2002, ApJ, 569, 1016
  * Feng et al. (2013) Feng, L., Wiegelmann, T., Su, Y., et al. 2013, ApJ, 765, 37
  * Georgoulis et al. (2012) Georgoulis, M. K., Titov, V. S., & Mikić, Z. 2012, ApJ, 761, 61
  * Glogowski et al. (2019) Glogowski, K., Bobra, M. G., Choudhary, N., Amezcua, A. B., & Mumford, S. J. 2019, Journal of Open Source Software, 4, 1614, doi: 10.21105/joss.01614
  * Guo et al. (2008) Guo, Y., Ding, M. D., Wiegelmann, T., & Li, H. 2008, ApJ, 679, 1629
  * Hagino & Sakurai (2004) Hagino, M., & Sakurai, T. 2004, PASJ, 56, 831
  * Hoeksema et al. (2014) Hoeksema, J. T., Liu, Y., Hayashi, K., et al. 2014, SoPh, 289, 3483
  * Kazachenko et al. (2017) Kazachenko, M. D., Lynch, B. J., Welsch, B. T., & Sun, X. 2017, ApJ, 845, 49
  * Kontogiannis et al. (2017) Kontogiannis, I., Georgoulis, M. K., Park, S.-H., & Guerra, J. A. 2017, SoPh, 292, 159, doi: 10.1007/s11207-017-1185-1
  * Kusano et al. (2012) Kusano, K., Bamba, Y., Yamamoto, T. T., et al. 2012, ApJ, 760, 31
  * Lee et al. (2018) Lee, E.-J., Park, S.-H., & Moon, Y.-J. 2018, SoPh, 293, 159, doi: 10.1007/s11207-018-1381-7
  * Leka & Barnes (2003a) Leka, K. D., & Barnes, G. 2003a, ApJ, 595, 1277
  * Leka & Barnes (2007) —. 2007, ApJ, 656, 1173
  * Li et al. (2022) Li, T., Sun, X., Hou, Y., et al. 2022, ApJL, 926, L14
  * Liu et al. (2017a) Liu, C., Deng, N., Wang, J. T. L., & Wang, H. 2017a, ApJ, 843, 104
  * Liu et al. (2013) Liu, Y. D., Luhmann, J. G., Lugaz, N., et al. 2013, ApJ, 769, 45
  * Liu et al. (2014a) Liu, Y. D., Richardson, J. D., Wang, C., & Luhmann, J. G. 2014a, ApJL, 788, L28
  * Liu et al. (2019) Liu, Y. D., Zhao, X., Hu, H., Vourlidas, A., & Zhu, B. 2019, ApJS, 241, 15
  * Liu et al. (2014b) Liu, Y. D., Luhmann, J. G., Kajdič, P., et al. 2014b, NatCo, 5, 3481
  * Lu et al. (2019) Lu, Z., Cao, W., Jin, G., et al. 2019, ApJ, 876, 133
  * Mason & Hoeksema (2010) Mason, J. P., & Hoeksema, J. T. 2010, ApJ, 723, 634
  * Moore et al. (2012) Moore, R. L., Falconer, D. A., & Sterling, A. C. 2012, ApJ, 750, 24
  * Moore et al. (1987) Moore, R. L., Hagyard, M. J., & Davis, J. M. 1987, SoPh, 113, 347
  * Petrie (2012a) Petrie, G. J. D. 2012a, ApJ, 759, 50
  * Petrie (2012b) —. 2012b, in Solar Dynamics and Magnetism from the Interior to the Atmosphere (Springer), 415–440
  * Ribeiro & Gradvohl (2021) Ribeiro, F., & Gradvohl, A. 2021, Astronomy and Computing, 35, 100468
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, SoPh, 275, 229
  * Schrijver et al. (2008) Schrijver, C., DeRosa, M., Metcalf, T., et al. 2008, ApJ, 675, 1637
  * Schrijver (2007) Schrijver, C. J. 2007, ApJ, 655, L117, doi: 10.1086/511857
  * Seehafer (1990) Seehafer, N. 1990, SoPh, 125, 219
  * Song & Zhang (2016) Song, Y. L., & Zhang, M. 2016, ApJL, 826, 173
  * Sun et al. (2012) Sun, X., Hoeksema, J. T., Liu, Y., et al. 2012, ApJ, 748, 77
  * Tarr et al. (2013) Tarr, L., Longcope, D., & Millhouse, M. 2013, ApJ, 770, 4
  * The SunPy Community et al. (2020) The SunPy Community, Barnes, W. T., Bobra, M. G., et al. 2020, ApJ, 890, 68, doi: 10.3847/1538-4357/ab4f7a
  * Török et al. (2014) Török, T., Leake, J. E., Titov, V. S., et al. 2014, ApJL, 782, L10, doi: 10.1088/2041-8205/782/1/l10
  * Vasantharaju et al. (2018) Vasantharaju, N., Vemareddy, P., Ravindra, B., & Doddamani, V. H. 2018, ApJ, 860, 58
  * Vemareddy & Wiegelmann (2014) Vemareddy, P., & Wiegelmann, T. 2014, ApJ, 792, 40
  * Wang et al. (2019) Wang, J., Liu, S., Ao, X., et al. 2019, ApJ, 884, 175, doi: 10.3847/1538-4357/ab441b
  * Wang et al. (2020) Wang, J., Zhang, Y., Webber, S. A. H., et al. 2020, ApJ, 892, 140, doi: 10.3847/1538-4357/ab7b6c
  * Wang et al. (2018) Wang, R., Liu, Y. D., Hoeksema, J. T., Zimovets, I. V., & Liu, Y. 2018, ApJ, 869, 90
  * Wang et al. (2022) Wang, R., Liu, Y. D., Yang, S., & Hu, H. 2022, ApJ, 925, 202
  * Wang et al. (2014) Wang, R., Liu, Y. D., Yang, Z., & Hu, H. 2014, ApJ, 791, 84
  * Wang et al. (2013) Wang, R., Yan, Y., & Tan, B. 2013, SoPh, 288, 507
  * Welsch et al. (2011) Welsch, B. T., Christe, S., & McTiernan, J. M. 2011, SoPh, 274, 131
  * Yang et al. (2013) Yang, X., Lin, G., Zhang, H., & Mao, X. 2013, ApJL, 774, L27
  * Yeates et al. (2008) Yeates, A. R., Mackay, D. H., & van Ballegooijen, A. A. 2008, ApJ, 680, L165
  * Yi et al. (2021) Yi, K., Moon, Y.-J., Lim, D., Park, E., & Lee, H. 2021, ApJ, 910, 8
  * Zhang (2001) Zhang, H. 2001, ApJ, 557, L71


# LHC signals of a BLSSM CP-even Higgs boson

A. Hammad1, S. Khalil1, S. Moretti2 1 Center for Fundamental Physics, Zewail
City of Science and Technology, 6 October City, Giza, Egypt.
2 School of Physics & Astronomy, University of Southampton, Highfield,
Southampton, UK.

###### Abstract

We study the scope of the Large Hadron Collider in accessing a neutral Higgs
boson of the $B-L$ Supersymmetric Standard Model. After assessing the
surviving parameter space configurations following the Run 1 data taking, we
investigate the possibilities of detecting this object during Run 2. For the
model configurations in which the mixing between such a state and the
discovered Standard Model-like Higgs boson is non-negligible, there exist
several channels enabling its discovery over a mass range spanning from
$\approx 140$ to $\approx$ 500 GeV. For a heavier Higgs state, with mass above
250 GeV (i.e., twice the mass of the Higgs state discovered in 2012), the
hallmark signature is its decay in two such 125 GeV scalars, $h^{\prime}\to
hh$, where $hh\to b\bar{b}\gamma\gamma$. For a lighter Higgs state, with mass
of order 140 GeV, three channels are accessible: $\gamma\gamma$, $Z\gamma$ and
$ZZ$, wherein the $Z$ boson decays leptonically. In all such cases,
significances above discovery can occur for already planned luminosities at
the CERN machine.

##  I Introduction

After the Higgs boson discovery at the Large Hadron Collider (LHC) during Run
1, a new era in particle physics has begun. While precision measurements of
the detected state as reported by the ATLAS and CMS collaborations (now also
including Run 2 data) confirm a Standard Model (SM)-like nature with a rather
light mass of $\approx 125$ GeV, significant effort is now being put in the
search for companion Higgs states, as any Beyond the SM (BSM) construct
embedding a Higgs mechanism is likely being non-minimal, i.e., it would
include new Higgs bosons in its spectrum. In the miyriad of BSM scenarios
available, a special place is held by models of Supersymmetry (SUSY), wherein
the lightest SM-like Higgs boson mass is naturally limited to be at the
Electro-Weak (EW) scale (say below $2M_{W}$) and where one also finds
additional (neutral) Higgs bosons. Thus, one may well be tempted to conclude
that a SUSY scenario may be behind the aforementioned data.

Amongst the many SUSY realisations studied so far, though, one really ought to
single out those that also offer explanations to other data pointing to BSM
physics, chiefly those indicating that neutrinos oscillate, hence that they
have mass. One is therefore well motivated in looking at the $B-L$
Supersymmetric Standard Model (BLSSM). The BLSSM is an extension of the time-
honoured Minimal Supersymmetric Standard Model (MSSM) obtained by adopting a
further $U(1)_{B-L}$ gauge group alongside the SM structure, i.e.,
$SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times U(1)_{B-L}$. (This requires an
additional Higgs singlet field to break the new $U(1)_{B-L}$ symmetry, in turn
releasing an additional $Z^{\prime}$ state as well.) The particle content of
the BLSSM, limited to its Higgs sector, includes three additional neutral
Higgs fields (henceforth $h^{\prime},H^{\prime}$ and $A^{\prime}$) with
respect to the MSSM ones (henceforth $h,H$ and $A$)111We conventionally denote
here by $A^{(^{\prime})}$ CP-odd Higgs states and by
$h^{(^{\prime})},H^{(^{\prime})}$ CP-even ones (the latter with $m_{h}<m_{H}$
and $m_{h^{\prime}}<m_{H^{\prime}}$). Notice that also two charge conjugated
states are present in both the MSSM and BLSSM, denoted by $H^{\pm}$..

The enriched Higgs sector of the BLSSM, with respect to the MSSM one, offers
the possibility of relieving the deadlock typical of the minimal SUSY model,
wherein a light SM-like Higgs state (the $h$ boson at $\approx 125$ GeV)
requires the other Higgs states ($H$ and $A$ in particular) to be much heavier
in comparison (and moderately coupled to SM matter fermions and gauge bosons).
This does not necessarily occurs in the BLSSM, as the $h^{\prime},H^{\prime}$
and $A^{\prime}$ states can have a singlet component sufficient to render them
very lightly mixed with the $h$ one, thereby allowing at the same time sizable
couplings to SM objects and the possibility of their mass, depending on the
Vacuum Expectation Value (VEV) of the Higgs singlet field, to be significantly
lighter than those of the MSSM-like $H$ and $A$ particles.

In fact, a natural configuration of the BLSSM is to find alongside the above
SM-like Higgs state another rather light physical Higgs boson, $h^{\prime}$,
also CP-even, with a mass $m_{h^{\prime}}\geq 135$ GeV. This fact was
exploited in Refs. Abdallah:2014fra ; Khalil:2015vpa ; Hammad:2015eca  to
explain potential Run 1 signals for another Higgs boson, i.e., $h^{\prime}$,
in the $h^{\prime}\to ZZ^{\ast}\to 4l$ (wherein a $2\sigma$ excess is
appreciable in the vicinities of 145 GeV Chatrchyan:2013mxa ),
$h^{\prime}\to\gamma\gamma$ (prompting a $2.9\sigma$ excess around 137 GeV
CERNAA ) and $h^{\prime}\to Z\gamma$ (yielding a $2\sigma$ excess around $140$
GeV Chatrchyan:2013vaa ) decay modes.

As new data are presently being collected at Run 2, we revisit here the scope
of the LHC in confirming or disproving the above hypothesis of additional
light Higgs boson signals. Furthermore, thanks to the higher energy and
luminosity afforded by the new CERN machine configuration, we also investigate
the possibilities of accessing a heavier $h^{\prime}$ state, with a mass up to
500 GeV or so. Our paper is organised as follows. In the next section, we
introduce the Higgs boson spectrum in the BLSSM. In Sects. III and IV we
describe our analysis of the light and heavy, respectively, mass range of the
$h^{\prime}$ state. We conclude in Sect. V.

##  II Higgs bosons in the BLSSM

The BLSSM model consists of, in addition to the MSSM particle content, two SM
singlet chiral Higgs superfields $\chi_{1,2}$ and three SM singlet chiral
superfields, $\nu_{i},i=1,2,3$ Khalil:2007dr . The Superpotential of this
model is given by

$W=Y_{u}\hat{Q}\hat{H}_{u}\hat{U}^{c}+Y_{d}\hat{Q}\hat{H}_{d}\hat{D}^{c}+Y_{e}\hat{L}\hat{H}_{d}\hat{E}^{c}+Y_{\nu}\hat{L}\hat{H}_{u}\hat{\nu}^{c}\\\
+Y_{\nu}\,\hat{\nu}\,\hat{l}\,\hat{H}_{u}.+\mu\hat{H}_{u}\hat{H}_{d}+\mu^{\prime}\hat{\chi}_{1}\hat{\chi}_{2}.$
(1)

The corresponding soft SUSY breaking terms and the details of the associated
spectrum can be found in Refs. Khalil:2007dr ; Staub:2015kfa . Note that the
$U(1)_{Y}$ and $U(1)_{B-L}$ gauge kinetic mixing can be absorbed in the
covariant derivative redefinition and, in this basis, one finds

$\displaystyle M_{Z}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{4}(g_{1}^{2}+g_{2}^{2})v^{2},$ (2) $\displaystyle
M_{Z^{\prime}}^{2}$ $\displaystyle=$ $\displaystyle g_{BL}^{2}v^{\prime
2}+\frac{1}{4}\tilde{g}^{2}v^{2},$ (3)

where $g_{BL}$ is the gauge coupling of $U(1)_{B-L}$ and $\tilde{g}$ is the
gauge coupling mixing between $U(1)_{Y}$ and $U(1)_{B-L}$. In addition,
$v=\sqrt{v^{2}_{1}+v^{2}_{2}}\simeq 246$ GeV, $v^{\prime}=\sqrt{v^{\prime
2}_{1}+v^{\prime 2}_{2}}\simeq{\cal O}(1)$ TeV are the VEVs of the Higgs
fields $H_{i}$ and $\chi_{i}$, respectively.

###  II.1 The spectrum

The neutral Higgs boson masses are obtained by making the usual redefinition
of the Higgs fields, i.e.,
$H_{1,2}^{0}={\frac{1}{\sqrt{2}}}(v_{1,2}+\sigma_{1,2}+i\phi_{1,2})$ and
$\chi_{1,2}^{0}={\frac{1}{\sqrt{2}}}(v^{\prime}_{1,2}+\sigma^{\prime}_{1,2}+i\phi^{\prime}_{1,2})$,
where $\sigma_{1,2}={\rm Re}H_{1,2}^{0}$, $\phi_{1,2}={\rm Im}H_{1,2}^{0}$,
$\sigma^{\prime}_{1,2}={\rm Re}\chi_{1,2}^{0}$ and $\phi^{\prime}_{1,2}={\rm
Im}\chi_{1,2}^{0}$. The real parts correspond to the CP-even Higgs bosons and
the imaginary parts correspond to the CP-odd Higgs bosons. Therefore, the
squared matrix of the BLSSM CP-even neutral Higgs fields at tree level, in the
basis $(\sigma_{1},\sigma_{2},\sigma^{\prime}_{1},\sigma^{\prime}_{2})$, is
given by

$M^{2}=\left(\begin{array}[]{cc}M^{2}_{h{H}}{}&~{}M^{2}_{hh^{\prime}}\\\ \\\
M^{2^{{T}}}_{hh^{\prime}}{}&~{}M^{2}_{h^{\prime}{H}^{\prime}}\end{array}\right),$
(4)

where $M^{2}_{h{H}}$ is the usual MSSM neutral CP-even Higgs mass matrix,
which leads to a SM-like Higgs boson with mass, at one loop level, of order
125 GeV and a heavy Higgs boson with mass $m_{H}\sim{\cal O}(1$ TeV). In
addition, the BLSSM matrix $M^{2}_{h^{\prime}{H}^{\prime}}$ is given by

$\displaystyle
M^{2}_{h^{\prime}{H}^{\prime}}=\left(\begin{array}[]{cc}m^{2}_{A^{\prime}}c^{2}_{\beta^{\prime}}+g^{2}_{BL}v^{\prime
2}_{1}&-\frac{1}{2}m^{2}_{A^{\prime}}s_{2\beta^{\prime}}-g^{2}_{BL}v^{\prime}_{1}v^{\prime}_{2}\\\
\\\
-\frac{1}{2}m^{2}_{A^{\prime}}s_{2\beta^{\prime}}-g^{2}_{BL}v^{\prime}_{1}v^{\prime}_{2}&m^{2}_{A^{\prime}}s^{2}_{\beta^{\prime}}+g^{2}_{BL}v^{\prime
2}_{2}\end{array}\right),$ (8)

where $c_{x}=\cos(x)$ and $s_{x}=\sin(x)$. Therefore, the eigenvalues of this
mass matrix are given by

$\displaystyle{m}^{2}_{h^{\prime},H^{\prime}}=\frac{1}{2}\Big{[}(m^{2}_{A^{\prime}}+M_{{Z^{\prime}}}^{2})\mp\sqrt{(m^{2}_{A^{\prime}}+M_{{Z^{\prime}}}^{2})^{2}-4m^{2}_{A^{\prime}}M_{{Z^{\prime}}}^{2}\cos^{2}2\beta^{\prime}}\;\Big{]}.$
(9)

If $\cos^{2}{{2}\beta^{\prime}}\ll 1$, one finds that the lightest $B-L$
neutral Higgs state is given by

$\displaystyle{m}_{h^{\prime}}\;{\simeq}\;\left(\frac{m^{2}_{A^{\prime}}M_{{Z^{\prime}}}^{2}\cos^{2}2\beta^{\prime}}{{m^{2}_{A^{\prime}}+M_{{Z^{\prime}}}^{2}}}\right)^{\frac{1}{2}}\simeq{\cal
O}(100~{}{\rm GeV}).$ (10)

The mixing matrix $M_{hh^{\prime}}^{2}$ is proportional to $\tilde{g}$ and,
for a gauge coupling $g_{BL}\sim|\tilde{g}|\sim{\cal O}(0.1)$, these off-
diagonal terms are about one order of magnitude smaller than the diagonal
ones. However, they are still crucial for generating interaction vertices
between the light BLSSM Higgs boson, $h^{\prime}$, and the MSSM-like Higgs
state, $h$.

The CP-even neutral Higgs mass matrix in Eq. (4) can be diagonalised by a
unitary transformation:

$\displaystyle{\Gamma}~{}M^{2}~{}\Gamma^{\dagger}={\rm
diag}\\{m_{h}^{2},m_{h^{\prime}}^{2},m_{H}^{2},m_{H^{\prime}}^{2}\\}.$ (11)

A numerical scan confirms that, while $\tan^{\prime}\beta\leq 1.2$, the
$h^{\prime}$ state can be the second Higgs boson mass whereas the other two
CP-even states $H,H^{\prime}$ are heavy. Also, the mixings $\Gamma_{ij}$ are
proportional to $\tilde{g}$ and they vanish (at tree level) if $\tilde{g}=0$.
In this regard, $h^{\prime}$ can be written in terms of gauge eigenstates as

$\displaystyle
h^{\prime}=\Gamma_{21}~{}\sigma_{1}+\Gamma_{22}~{}\sigma_{2}+\Gamma_{23}~{}\sigma^{\prime}_{1}+\Gamma_{24}~{}\sigma^{\prime}_{2}.$
(12)

Thus, the couplings of the $h^{\prime}$ with up- and down-quarks are given by

$g_{h^{\prime}u\bar{u}}=-i\frac{m_{u}\times\Gamma_{22}}{\upsilon\sin\beta},$
(13)
$g_{h^{\prime}d\bar{d}}=-i\frac{m_{d}\times\Gamma_{21}}{\upsilon\cos\beta}.$
(14)

Similarly, one can derive the $h^{\prime}$ couplings with the ${W^{+}}{W^{-}}$
and $Z{Z}$ gauge boson pairs:

$\displaystyle g_{{}_{h^{\prime}WW}}$ $\displaystyle=$
$\displaystyle{i}~{}g_{{2}}M_{W}\left(\Gamma_{{22}}\sin\beta+\Gamma_{{21}}\cos\beta\right),$
$\displaystyle g_{{}_{h^{\prime}ZZ}}$ $\displaystyle=$
$\displaystyle\frac{i}{2}\Big{[}4g_{BL}\sin^{2}{\theta^{\prime}}\left(v^{\prime}_{1}\Gamma_{22}+v^{\prime}_{2}\Gamma_{21}\right)+\left(v_{2}\Gamma_{22}+v_{1}\Gamma_{21}\right)\left(g_{z}\cos{\theta^{\prime}}-\tilde{g}\sin{\theta^{\prime}}\right)^{2}\Big{]},~{}~{}$
(15)

where $g_{z}=\sqrt{g_{1}^{2}+g_{2}^{2}}$ and $\theta^{\prime}$ is the mixing
angle between $Z$ and $Z^{\prime}$. Since $\sin{\theta^{\prime}}\ll 1$ (as per
experimental constraints), the coupling of the $h^{\prime}$ with $ZZ$,
$g_{{}_{h^{\prime}ZZ}}$, will be as follows:

$\displaystyle g_{{}_{h^{\prime}ZZ}}\simeq
i~{}g_{z}\,M_{Z}\left(\Gamma_{22}\sin{\beta}+\Gamma_{21}\cos{\beta}\right).$
(16)

Figure 1: Decomposition of the BLSSM Higgs boson, $h^{\prime}$, and the SM-
like Higgs, $h$, versus $M_{h^{\prime}}$.

In Fig. 1 we show the $h^{\prime}$ (in $\Gamma_{21},\Gamma_{22}$) and $h$ (in
$\Gamma_{11},\Gamma_{12}$) decompositions. Note that, if $\tilde{g}=0$, the
coupling of the BLSSM lightest Higgs boson with the SM particles vanishes at
tree level and is very suppressed ($\sim{\cal{O}}(10^{-6})$) at loop level.
Here we choose a parameter space such that the lightest chargino is rather
light, $M_{\chi^{\pm}}=120$ GeV, so as to enhance the SUSY contributions to
the Higgs decays into $\gamma\gamma$ and $Z\gamma$, namely, we consider a low
$\tan\beta$ between $1.1$ and $5$ and $\mu$ and $M_{2}$ between $100$ and
$300$ GeV, while other SUSY mass and trilinear parameters are assumed to be of
order few TeV. It is worth mentioning that the dominant decomposition for the
SM-like Higgs state is $\Gamma_{12}\sim{\cal O}(1)$, which is equivalent to
$\sin\beta\sim{\cal O}(1)$ in the MSSM, and that the light BLSSM Higgs,
$h^{\prime}$, is dominated by $\Gamma_{23}$ and $\Gamma_{24}\sim{\cal
O}(0.5)$.

We display in Fig. 2 the Branching Ratios (BRs) of $h^{\prime}$ into all its
possible decay channels, for non-zero $\tilde{g}$, including $gg$,
$\gamma\gamma$ and $Z\gamma$ that are induced at one loop level. A few remarks
on this figure are in order: $(i)$ for $m_{h^{\prime}}\geq 200$ GeV,
$h^{\prime}$ decays are dominated by the $W^{+}W^{-}$ and $hh$ channels;
$(ii)$ in the BLSSM the BR$(h^{\prime}\to Z\gamma)$ is typically larger than
the BR$(h^{\prime}\to\gamma\gamma)$, unlike the MSSM and SM where it is the
other way around.

Figure 2: The BRs of $h^{\prime}$ versus $M_{h^{\prime}}$ for
$0.1\leq\tilde{g}\leq 0.25$ and $g_{BL}=0.5$.

###  II.2 Implementation and simulation

The Higgs production modes included in our forthcoming numerical analysis are
gluon-gluon Fusion (ggF), which induce around $90\%$ of the total cross
section (hereafter denoted by $\sigma$), while Vector-Boson Fusion (VBF),
Higgs-strahlung (VH) and associated production with top-quarks (ttH)
contribute with around $10\%$. The data analyses in these channels are based
on an integrated luminosity of $20$ fb-1 at $\sqrt{s}=7,8$ TeV and expected to
rely upon from $100$ to ${\cal O}(1000)$ fb-1 at $\sqrt{s}=13$ TeV. The
magnitude of the signal is usually expressed via the “signal strength”
parameters, defined as

$\displaystyle\mu_{{}_{XY}}$ $\displaystyle=$ $\displaystyle\frac{\sigma(pp\to
h^{\prime}\to XY)}{\sigma(pp\to h\to XY)^{\rm SM}}=\frac{\sigma(pp\to
h^{\prime})}{\sigma(pp\to h)^{\rm SM}}\times\frac{{\rm BR}(h^{\prime}\to
XY)}{{\rm BR}(h\to XY)^{\rm SM}}.$ (17)

Herein, the $h^{\prime}$ in the numerator is indeed the lightest BLSSM CP-even
state and the $h$ in the denominator is the SM Higgs boson with mass 125 GeV
whereas in both cases the cross section is intended as computed inclusively
(i.e., over the ggF, VBF, VH and ttH modes222Hereafter, the bulk of the
production rates will be due to the first two channels.).

For the implementation of the BLSSM we used SARAH Staub:2013tta  and SPheno
Porod:2011nf  to build the model. For loop induced channels we linked it with
CP-SuperH Lee:2012wa . The matrix-element calculation and events generation
were derived by MadGraph Alwall:2011uj . We then used Pythia Sjostrand:2007gs
to simulate the initial and final state radiation, fragmentation and
hadronisation effects. For detector simulation we passed the Pythia output to
Delphes deFavereau:2013fsa . For data analysis, we used MadAnalysis5
Conte:2012fm .

In our scans, for the computation of the signal strength distributions in the
next section, we consider the following regions of the parameter space:

$\displaystyle m_{0}$ $\displaystyle=$ $\displaystyle 1-3~{}{\rm
TeV},M_{3}=3~{}{\rm TeV},M_{2}=120-300~{}{\rm GeV},M_{1}=100-500~{}{\rm
GeV},\tan\beta=5,$ $\displaystyle\tan\beta^{\prime}$ $\displaystyle=$
$\displaystyle 1.15,|A_{0}|=1.5-3~{}{\rm TeV},\mu=100-350~{}{\rm
GeV},|\tilde{g}|=0.1-0.25,g_{BL}=0.5.$ (18)

In addition, in the upcoming event generation analyses, the following
benchmark point is assumed:

$\displaystyle m_{\chi_{1}^{+}}=120~{}{\rm GeV},\mu=120~{}{\rm
GeV},\tan\beta=5,\tan\beta^{\prime}=1.15,\tilde{g}=-0.24,g_{BL}=0.5,$ (19)

while all other SUSY particles are of order TeV. This benchmark point is
consistent with current theoretical and experimental limits, as we determined
through an independent program checked against specialised literature. It is
worth pointing out that light $\mu$ and chargino mass are crucial for
enhancing the SUSY contributions to $h\to\gamma\gamma$ and $h\to Z\gamma$
simultaneously. Finally, notice that the $h^{\prime}$ masses considered below
(140, 300, 350 and 480 GeV) are all accessible through the inputs in Eq. (19),
upon suitable adjustments of the Higgs potential parameters.

##  III Search for a heavy BLSSM Higgs boson at the LHC

In this section we analyse possible signatures of the lightest genuine BLSSM
scalar boson $h^{\prime}$ when it is rather heavy, with mass between 300 GeV
and 1 TeV, at Run 2 of the LHC. Fig. 2 shows that the decay channels available
to the $h^{\prime}$ state are the same as those of the SM-like $h$ one, with
the notable exception of the former decaying into (pairs of) the latter, i.e.,
$h^{\prime}\to hh$. The corresponding BR can be in fact the dominant one, once
its threshold is open. It is therefore the distinctive feature of a heavy
$h^{\prime}$ whenever $m_{h^{\prime}}\geq 2m_{h}$.

ATLAS ATLAS:2014rxa  and CMS CMS:2014eda  have both recently searched for $hh$
signals decaying to a $4b$ final state. However, it turned out to be a
significant challenge to distinguish the emerging signature, made of of four
$b$-jets in the final state, from the huge multi-jet QCD background. In fact,
the sensitivity achieved by the LHC experiments was rather poor and results
obtained were consistent with the SM. We shall nonetheless attempt extracting
this signal, so as to compare the scope of detecting it at Run 2 versus what
has been achieved at Run 1.

The decay $h^{\prime}\to hh\to\gamma\gamma b\bar{b}$, which has been
experimentally analysed in Refs. Aad:2014yja ; CMS:2014ipa , may prove to be
the best way to probe a heavy $h^{\prime}$ of the BLSSM, since the problem of
a suppressed $h\to\gamma\gamma$ decay is offset by the fact that both
$h^{\prime}\to hh$ and $h\to b\bar{b}$ are the dominant decays of the two
Higgs states concerned. The aforementioned searches were performed on the
$\sqrt{s}=8$ TeV data set corresponding to an integrated luminosity of
$\approx 20~{}{\rm fb}^{-1}$. Following these, the ATLAS collaboration
observed five excess events (above and beyond the expected SM yield) within a
mass windows from 260 to 500 GeV, which represent an excess of $2.4\sigma$,
with an intriguing $p_{0}$-value (local probability of compatibility with the
background) $\sim 10^{-3}$ at 300 GeV, which corresponds to $3.0~{}\sigma$
Aad:2014yja . In contrast, CMS reported that searches within the mass region
from 260 GeV to 1100 GeV were consistent with expectations from SM processes
CMS:2014ipa . Needless to say then, we will thoroughly investigate this
signature too at the upcoming Run 2.

Before proceeding to doing so in two separate subsections, let us start by
explaining how such large decay rates for $h^{\prime}\to hh$ can occur in the
BLSSM. Herein, the scalar trilinear coupling between $h^{\prime}$ and $hh$ is
given by

$\displaystyle\lambda_{h^{\prime}hh}^{\rm
BLSSM}=\frac{-i\tilde{g}g_{{}_{BL}}}{4}\Gamma_{i2}^{2}\left(2v^{\prime}_{2}\Gamma_{24}-v^{\prime}_{1}\Gamma_{23}\right).$
(20)

Here we have assumed, as advocated in the previous section, that
$\Gamma_{12}\gg\Gamma_{11,13,14}$ and $\Gamma_{23,24}\gg\Gamma_{21,22}$. This
should be compared with the MSSM trilinear scalar coupling

$\displaystyle\lambda_{Hhh}^{\rm
MSSM}=-i\frac{g_{1}^{2}+g_{2}^{2}}{4}v\left[2\sin
2\alpha\sin(\beta+\alpha)-\cos 2\alpha\cos(\beta+\alpha)\right],$ (21)

for which, when $\sin\beta>\cos\beta$ and assuming the decoupling limit where
$\alpha\sim\beta$, one finds

$\lambda_{Hhh}^{\rm MSSM}=-i\frac{g_{1}^{2}+g_{2}^{2}}{4}v\sin^{3}\beta.$ (22)

Also note that the ${Hhh}$ coupling is modified in the BLSSM with respect to
the MSSM and takes the form

$\displaystyle\lambda_{Hhh}^{\rm
BLSSM}=\frac{i}{4}(g_{1}^{2}+\tilde{g}^{2}+g_{2}^{2})\Gamma_{31}\left(v_{d}\Gamma_{12}^{2}+2v_{u}\Gamma_{12}\Gamma_{11})\right).$
(23)

It is clear that $\lambda_{h^{\prime}hh}^{\rm BLSSM}\propto
v^{\prime}_{1,2}\sim{\cal O}(1)$ TeV is much larger than the coupling ${Hhh}$
in either SUSY model, which is of order of the EW scale. Therefore, one would
expect that the decay rate of $h^{\prime}$ $\to$ $hh$ is always much larger
than that of $H\to hh$.

###  III.1 The $hh\to 4b$ decays of a heavy BLSSM Higgs boson

The total cross section for the aforementioned $4b$ final state is given by

$\displaystyle\sigma(pp\to h^{\prime}\to hh\to 4b)=\sigma(pp\to
h^{\prime})\times{\rm{BR}}(h^{\prime}\to hh)\times{\rm{BR}}(h\to
b\bar{b})^{2},$ (24)

and is dominated by ggF which is in turn obtained as (for a CM energy of 13
TeV)

$\sigma(pp\to h)\times\Gamma_{22}^{2}\simeq{\cal O}(1)~{}{\rm pb}$ (25)

while, for $m_{h^{\prime}}\simeq 350$ GeV, the BR$(h^{\prime}\to hh)\sim 0.5$
and the BR$(h\to b\bar{b})\sim 0.6$, as can be seen from Fig. 2. Thus, one
finds that $\sigma(pp\to h^{\prime}\to hh\to 4b)$ in the BLSSM $\sim
10^{-1}~{}{\rm pb}$. Altough the high total cross section, the huge
contribution from background b-jet radiation exceed the signal, so that the
associated events would not appear as significant over the SM background. This
conclusion is confirmed by Fig. 3, where we show the number of events of
signal with its irreducible background as a function of the invariant mass of
the four $b$-jets, $M_{4b}$. Note that we used the $b$-tagging algorithm
included in MadAnalysis Conte:2012fm , so that a jet is identified as
originating from a $b$-quark when it can be matched to it once it lies within
a cone of radius certain $R$ around one of the parton-level $b$-quarks, this
yielding an efficiency of about $65\%$.

Here, we considered the cuts applied in Aad:2015uka : i.e., candidate events
are required to have at least four $b$-tagged jets, each with $p_{T}\geq 40$
GeV and separated by a cone of $\Delta R=1.5$. However, as can be seen from
the plot, the signal is well below the background. The highest background
contribution comes from a muti-jets final state, followed by $t\bar{t}$
production and (semi-)hadronic (anti)top decays which gives about $22\%$ of
the noise while the reducible background contributions come from ‘$Z$ \+
jets’, $ZZ$ and $Zh$ and are found to contribute less than $1\%$.

The signal distribution is presented for $m_{h}^{\prime}\approx 2m_{t}\approx
350$ GeV, which is in fact the worse case scenario, as this is where the
$t\bar{t}$ background peaks in $M_{4b}$. However, we have tried different
$m_{h^{\prime}}$ values, to no avail, in the mass range from 300 GeV to 1 TeV.
The signal would never be accessible, neither with standard nor with upgraded
luminosities.

Figure 3: Number of signal events for $h^{\prime}\to hh\to 4b$ decays (red)
induced by ggF and VBF versus the $4b$ invariant mass at $\sqrt{s}=13$ TeV
after 100 fb-1 of luminosity alongside the $t\bar{t}$ background (blue). (The
huge multi-jet background, which is given in Ref. Aad:2015uka , is not shown.)
Here, $m_{h^{\prime}}=350$ GeV.

###  III.2 The $hh\to b\bar{b}\gamma\gamma$ decays of a heavy BLSSM Higgs
boson

Now we turn to the process $pp\to h^{\prime}\to hh\to\gamma\gamma b\bar{b}$.
Although this mode has smaller cross section than $\sigma(pp\to h^{\prime}\to
hh\to 4b)$, it is more promising due to the clean di-photons trigger with
excellent mass resolution and low background contamination. This is confirmed
in Fig. 4, where the number of signal events is displayed versus the
background as a function of the invariant mass $M_{\gamma\gamma bb}$ for two
examples of $h^{\prime}$ masses: $m_{h}^{\prime}=300$ GeV and
$m_{h^{\prime}}=480$ GeV.

Figure 4: Number of signal events for $h^{\prime}\to\gamma\gamma b\bar{b}$
decays (red) induced by ggF and VBF versus the $\gamma\gamma b\bar{b}$
invariant mass at $\sqrt{s}=13$ TeV after 100 fb-1 of luminosity alongside the
two dominant $\gamma\gamma$ (blue) and $Zh$ (green) backgrounds. Their sum is
also shown as data points. Here, $m_{h}^{\prime}=300$ and $480$ GeV.

The background to this process can be classified into two categories:
background events containing a real Higgs boson decay, $h\to\gamma\gamma$ and
$h\to bb$, and the continuum background of events not containing a Higgs
boson. The continuum contribution in the signal region is split between events
with two photons and events with a single photon in association with a jet
faking the second photon. Further, the two $b$-tagged jets include real heavy-
flavour jets as well as mis-tagged light-flavour jets and gluons. The
contribution from di-leptonic decays of $t\bar{t}$ events where two electrons
fake the two photons is roughly 10$\%$ of the total background. The
contribution from other processes, like leptonic decays of di-gauge bosons
where two electrons fake the two photons and the Higgs boson comes associated
with a $W/Z$, is negligible. In our analysis, we adopt the following
acceptance cuts in transverse momentum, pseudorapidity of and separation
amongst the photons and jets:

  1. 1.

the pseudorapidity $\eta$ of the two photons must fall within the geometric
acceptance of the detector for photons, $|\eta|\leq 2.4$;

  2. 2.

the ratio between the transverse momenta of the leading and subleading photon
must be $\geq 0.25$;

  3. 3.

jets are required to fall within the tracker acceptance of $|\eta|\leq 2.5$
with transverse momentum $p_{T}\geq 35$ GeV.

After our preselection is enforced, already at standard luminosity of Run 2,
the signal is clearly visible above all backgrounds, both at 300 and 480 GeV,
thereby enabling one to declare discovery of a Higgs-to-two-Higgs signal as
well as circumstantial evidence of a BLSSM decay chain of the type
$h^{\prime}\to hh$. In order to eventually profile the latter though, the
simultaneous reconstruction of the two $h$ resonances and of the $h^{\prime}$
one is a pre-requisite. To this end, in Fig. 5, we also show the mass
reconstruction of the two SM-like Higgs boson masses, in the two channels
$h\to\gamma\gamma$ and $b\bar{b}$, against the backdrop of the SM noise. From
the corresponding distributions, a clear element emerges that characterises
this signature is very promising, i.e., the very efficient reconstriction of
$m_{h}\approx 125$ GeV from the di-photon pair, from which is evident the
strong background suppression which can be achieved. In contrast, this is not
true in the case of $b\bar{b}$ decays, as here the background remains
overwhelming above the signal (implicitly also explaining the reduced
sensitivity of the fully hadronic $4b$ signal previously considered, where jet
combinatorics would further play a significant role in degrading the quality
of it). Furthermore, notice that the quality of the mass reconstruction is not
dramatically different for $m_{h^{\prime}}=300$ and $480$ GeV.

Figure 5: Number of signal events for $h^{\prime}\to\gamma\gamma b\bar{b}$
decays (red) induced by ggF and VBF versus the $\gamma\gamma$ (left) and
$b\bar{b}$ (right) invariant mass at $\sqrt{s}=13$ TeV after 100 fb-1 of
luminosity alongside the total background (blue). Their sum is also shown as
data points. Here, $m_{h}^{\prime}=300$ and $480$ GeV. Only the acceptance
cuts described in the text are used here.

In the light of the mass distributions just discussed, one can attempt a more
refined signal selection against the continuum noise. In Tab. 2 we show the
number of events for signal and continuum background after each cut mentioned
therein and Fig. 6 shows the final number of events versus the background
after all cuts are applied. It is clear from this plot that the final result
is an almost background-free $M_{\gamma\gamma b\bar{b}}$ distribution neatly
pointing to the value of the $h^{\prime}$ mass, for values between 300 and 480
GeV. It is not surprising then, in the end, significances for the signal can
be extremely large, as seen in Fig. 7, for any $m_{h^{\prime}}$ value, after a
final sampling in $M_{\gamma\gamma b\bar{b}}$ is exploited. Notice that, here,
both reducible and irreducible backgrounds are accounted for in the
calculation.

Applied cut | Signal, $m_{h^{\prime}}=300$ | Signal, $m_{h^{\prime}}=480$ | Continuum background
---|---|---|---
After acceptance cuts | 626 | 237 | 4758
$M_{\gamma\gamma}\leq 135$ GeV | 625 | 234 | 4375
$M_{\gamma\gamma}\geq 115$ GeV | 616 | 223 | 182
$M_{b\bar{b}}\leq 145$ GeV | 536 | 210 | 98
$M_{b\bar{b}}\geq 105$ GeV | 351 | 86 | 30
Table 1: Signal (for two $h^{\prime}$ mass values) and continuum background
events in the $\gamma\gamma b\bar{b}$ channel as a function of several mass
selection cuts. The energy is $\sqrt{s}=13$ TeV whereas the luminosity is 100
fb-1. Figure 6: Number of signal events for $h^{\prime}\to\gamma\gamma
b\bar{b}$ decays (red) induced by ggF and VBF versus the $\gamma\gamma$ (left)
and $b\bar{b}$ (right) invariant mass at $\sqrt{s}=13$ TeV after 100 fb-1 of
luminosity alongside the total background (blue). Here, $m_{h}^{\prime}=300$
and $480$ GeV. Also the selection cuts of Tab. 2 are used here. Figure 7:
Left: Significance of the $h^{\prime}\to\gamma\gamma b\bar{b}$ signal (for
$m_{h^{\prime}}=300$ and 480 GeV) versus the luminosity (black). Right: Number
of events for signal and background for variable luminosity (red). Data are
produced at $\sqrt{s}=13$ TeV and the points correspond to an integrated
luminosity of 100, 300, 1000 and 3000 fb-1. Notice that event rates are
computed after the acceptance cuts described in the text and the mass
selections of Tab. 2. The $M_{\gamma\gamma b\bar{b}}$ mass windows used for
the calculation is 50 GeV for $m_{h^{\prime}}=300$ GeV and 100 GeV for
$m_{h^{\prime}}=480$ GeV.

##  IV Search for a light BLSSM Higgs boson at the LHC

In this section we briefly revisit the possible signatures of a light BLSSM
Higgs boson $h^{\prime}$ (with mass $m_{h^{\prime}}\approx 140$ GeV) at the
LHC. As emphasised in Refs. Abdallah:2014fra ; Hammad:2015eca ; Khalil:2015vpa
, this particle can be probed in one of the following channels:
$\gamma\gamma$, $Z\gamma$ and $ZZ$. We review these in the three upcoming
subsections.

###  IV.1 The $\gamma\gamma$ decays of a light BLSSM Higgs boson

The coupling of a Higgs boson with di-photons is induced by loops of charged
particles. In the SM, these loops are mediated by the $W$ gauge boson and top-
quark. In SUSY models, the $h\gamma\gamma$ triangle coupling contains
additional loops of charged particles: charged Higgses $H^{\pm}$ , squarks
$\tilde{q}$, sleptons $\tilde{\ell}^{\pm}$ and charginos $\chi^{\pm}$. Since
the Higgs boson coupling with SUSY particles are not proportional to their
masses their contributions decouple for high masses. In this paper, we focus
on the cases of light charginos, $\chi^{\pm}_{1}$, enhancements, since they
can increase the $h\gamma\gamma$ amplitude squared up to $30\%$
Belyaev:2013rza ; Hemeda:2013hha  (i.e., the sfermions and charged Higgs
bosons are assumed to be heavy).

The Higgs decay into di-photons provides a clean final-state topology which
allows for the mass to be reconstructed with high precision. The partial decay
width of the lightest BLSSM Higgs boson into di-photons is given by

$\Gamma(h^{\prime}\to\gamma\gamma)=\frac{G_{\mu}\alpha^{2}m^{3}_{h^{\prime}}}{128\sqrt{2}\pi^{3}}\left|{{A}}_{t}+{{A}_{W}}+A_{H^{\pm}}+A_{\tilde{f}}+{{A}_{\chi^{\pm}}}\right|^{2},$
(26)

where the amplitudes $A_{f,W,H^{\pm},{\tilde{f}},\chi^{\pm}}$ can be found in
Djouadi:2005gj . In Fig. 8 we show the signal strength of $gg\to
h^{\prime}\to\gamma\gamma$ for 110 GeV $<m_{h^{\prime}}<$ $150$ GeV. We also
include the di-photon signal strengths of the SM-like Higgs, $h$, in the MSSM
and BLSSM, in addition to the MSSM-like heavy Higgs, $H$. It is interesting to
note that the BLSSM results for both $h$ and $h^{\prime}$ are matching the
observed data at Run 1, whereas the signal strength of the heavy Higgs in the
MSSM, $H$, is quite suppressed and cannot easily account for these
observations.

Figure 8: Signal strength of the lightest and next-to-lightest Higgs bosons in
the BLSSM (in blue and red, respectively) in the $\gamma\gamma$ channel. For
comparison, we also include the signal strength of the lightest and next-to-
lightest Higgs bosons in the MSSM (in cyan and black, respectively). The $1$
and $2\sigma$ confidence intervals are extracted from data collected during
Run 1 with the observed exclusion limit as given in Aad:2014eha  is also
included.

The number of events for $h^{\prime}\to\gamma\gamma$ as function of the di-
photon invariant mass is presented in Fig. 9, for a Center-of-Mass (CM) energy
$\sqrt{s}=13\ {\rm TeV}$ and integrated luminosity = $100~{}{\rm fb}^{-1}$.
Here we choose the input parameters such that the SM-like Higgs boson has a
mass $m_{h}=125\ {\rm GeV}$ and the lightest genuinely BLSSM Higgs state has a
mass $m_{h^{\prime}}\sim 140\ {\rm GeV}$. The dominant backgrounds consist of
an irreducible fraction from prompt di-photon production and a reducible one
from $\gamma+$ jet and di-jet events where one or more of the objects
reconstructed as a photon corresponds to a jet, according to CMS “fake rates”.
It is also worth mentioning that here we consider all cuts applied in the CMS
analysis of Ref. Aad:2014eha : i.e, the photon candidates are collected within
$|\eta^{\gamma}|\leq 2.5$ with transverse momentum $p_{T}^{\gamma}\geq 20$
GeV. The production is considered here as induced from both ggF and VBF (as at
higher energies the latter mode grows in importace relatively to the former)
and yield both a $h$ and $h^{\prime}$ state. As can be seen from this figure,
the peak at $\sim 140$ GeV is greatly overwhelmed by the background after 100
fb-1, yet accessible with additional luminosity, as shown in Fig. 10.

Figure 9: Number of signal events for $h$ and $h^{\prime}\to\gamma\gamma$
decays (red and green, respectively) induced by ggF and VBF versus the
$\gamma\gamma$ invariant mass at $\sqrt{s}=13$ TeV after 100 fb-1 of
luminosity alongside the total background (blue). Their sum is also shown as
data points. Figure 10: Left: Significance of the $h^{\prime}\to\gamma\gamma$
signal (for $m_{h^{\prime}}=140$ GeV) versus the luminosity (black). Right:
Number of events for signal and background for variable luminosity (red). Data
are produced at $\sqrt{s}=13$ TeV and the points correspond to an integrated
luminosity of 100, 300, 1000 and 3000 fb-1. Notice that event rates are
computed after the cut $|m_{\gamma\gamma}-m_{h^{\prime}}|<10$ GeV.

###  IV.2 The $Z(\to\ell^{+}\ell^{-})\gamma$ decays of a light BLSSM Higgs
boson

Despite its small BR, the LHC experiments are currently sensitive to this
channel and will be so more and more as luminosity accrues. Precisely because
the SM rate in this decay channel is small, ATLAS and CMS may access BSM
physics through it, owing to the fact that the partial width can increase
sizeably in presence of additional loops of charged particles, just like in
the $h^{\prime}\to\gamma\gamma$ channel. The partial decay width of the
lightest BLSSM Higgs boson into $Z\gamma$ is given by

$\Gamma(h^{\prime}\to
Z\gamma)=\frac{G^{2}_{f}\alpha^{2}M^{2}_{W}M^{3}_{h^{\prime}}}{64\pi^{4}}\left(1-\frac{M^{2}_{Z}}{M^{2}_{h^{\prime}}}\right)^{3}\left|A_{f}+A_{W}+A_{H^{\pm}}+A_{\tilde{f}}+A_{\chi^{\pm}}\right|^{2},$
(27)

where the amplitudes $A_{f,W,H^{\pm},{\tilde{f}},\chi^{\pm}}$ can be found in
Djouadi:1996yq . As discussed in Hammad:2015eca , due to the mixing in the
sfermion and chargino sectors, the diagonal coupling only enhances the
$h^{\prime}\to\gamma\gamma$ channel, while the fact that the $Z$ boson has
both vector and axial vector quantum numbers makes both diagonal and off-
diagonal couplings of sfermions and charginos contribute to the $h^{\prime}\to
Z\gamma$ channel. As in $h^{\prime}\to\gamma\gamma$, we focus here on a light
chargino in order to enhance the $h^{\prime}\to Z\gamma$ amplitude.

In Fig. 11 we show that the signal strength of the $h^{\prime}$, $h$ (both in
the MSSM and BLSSM) and $H$ decays to $Z\gamma$ for $m_{h^{\prime},H}$ around
$140$ GeV (as usual, $m_{h}=125$ GeV), with the 1 and $2\sigma$ confidence
intervals extracted from data collected during Run 1 and with the observed
exclusion limit as given in Aad:2014fia . As can be seen, again, the BLSSM
results for both $h$ and $h^{\prime}$ match with the observed data rather well
whereas the signal strength of the heavy Higgs in the MSSM, $H$, as expected,
is quite suppressed, hence unable to reach out to current experimental
results.

Figure 11: Signal strength of the lightest and next-to-lightest Higgs bosons
BLSSM (in blue and red, respectively) in the $Z\gamma$ channel. The signal
strength of lightest and next-to-lightest Higgs bosons in the MSSM are given
in cyan and black points, respectively. The signal strength of the lightest
and next-to-lightest Higgs bosons in the MSSM are given in cyan and black,
respectively. The $1$ and $2\sigma$ confidence intervals are extracted from
data collected during Run 1 with the observed exclusion limit as given in
Aad:2014fia  is also included.

The distribution of the ‘di-lepton + photon’ (we assume
$Z\to\ell^{+}\ell^{-}$, $\ell=e,\mu$) invariant mass is presented in Fig. 12
for the signal and background, where the dominant components of the latter
consist of the irreducible contribution from $Z\gamma$ production, the
reducible one from final state radiation in the neutral Drell-Yan process and
‘$Z$ \+ jets’ processes where a jet is misidentified as a photon. Here the
cuts applied are as in Ref. Aad:2014fia , i.e.:

  1. 1.

the photon pseudorapidity must be $|\eta^{\gamma}|\leq 2.5$;

  2. 2.

the photon transverse momentum must be $p_{T}^{\gamma}\geq\ 25$ GeV;

  3. 3.

the di-lepton invariant mass must be $85$ GeV $\leq M_{\ell^{+}\ell^{-}}\leq
95$ GeV;

  4. 4.

the ‘di-lepton + photon’ invariant mass must be $130$ GeV $\leq
M_{\ell^{+}\ell^{-}\gamma}\leq 150$ GeV.

The cut flow results are found in Tab. 2. The selection (based on Run 1 cuts)
remains effective at Run 2 as well, since already at standard luminosity there
could already be an evidence of the $h^{\prime}\to Z\gamma$ signal in the
BLSSM. The line-shape of the signal, initially swamped by the background (see
left-hand side of Fig. 12), would also be very distinctive after the selection
is enforced (see right-hand side of Fig. 12). As the luminosity at Run 2
accumulates, the evidence will eventually turn into clear discovery (see Fig.
13).

| Signal ($S$) | Background ($B$) | $\frac{S}{\sqrt{S+B}}$
---|---|---|---
Before cuts | 200 | 18828 | 1.44
$p_{T}^{\gamma}\geq 25$ GeV | 180 | 6490 | 2.2
$85$ GeV $\leq M_{\ell^{+}\ell^{-}}\leq 95$ GeV | 172 | 4500 | 2.5
$130$ GeV $\leq M_{\ell^{+}\ell^{-}\gamma}\leq 150$ GeV | 170 | 3822 | 2.7

Table 2: Signal and background events in the $Z\gamma$ channel assuming
electron and muon decays of the $Z$ boson as a function of the selection cuts
detailed in the text. The energy is $\sqrt{s}=13$ TeV whereas the luminosity
is 100 fb-1.

Figure 12: Number of signal events for $h$ and $h^{\prime}\to
Z(\to\ell^{+}\ell^{-})\gamma$ decays ($\ell=e,\mu$) (red)) induced by ggF and
VBF versus the $\ell^{+}\ell^{-}\gamma$ invariant mass at $\sqrt{s}=13$ TeV
after 100 fb-1 of luminosity alongside the total background (yellow). Their
sum is also shown as data points. Left(Right): Before(After) the cuts in the
text are applied. Figure 13: Left: Significance of the $h^{\prime}\to
Z(\to\ell^{+}\ell^{-})\gamma$ signal (for $m_{h^{\prime}}=140$ GeV and
$\ell=e,\mu$) versus the luminosity (black). Right: Number of events for
signal and background for variable luminosity (red). Data are produced at
$\sqrt{s}=13$ TeV and the points correspond to an integrated luminosity of
100, 300, 1000 and 3000 fb-1. Notice that event rates are computed after the
cuts described in the text.

###  IV.3 The $ZZ(\to 4\ell)$ decays of a light BLSSM Higgs boson

The four leptons final state through the Higgs decay via pairs of $Z$ bosons
is the most significant channel for Higgs detection, yet it may not be the
most sensitive one to BSM effects, as its leading contribution occurs at tree
level, so that mixing effects of the SM-like boson with additional Higgs boson
states typically drive the BSM deviations. It was however one of the channels
where an anomaly at around 140 GeV appeared following the Run 1 analyses, as
intimated. In the MSSM, as mentioned above, in order to keep the signal
strength of the lightest Higgs boson $h$ consistent with the observed data,
one is constrained to the decoupling region, where at $M_{A}\gg M_{Z}$ and the
Higgs mixing angle $\alpha\sim\beta-\frac{\pi}{2}$. Therefore, the coupling of
the heaviest MSSM CP-even Higgs boson, $H$, with the SM gauge bosons is very
suppressed. In the case of the BLSSM, $\tilde{g}$ plays an important role in
enhancing both the first and the second lightest CP-even Higgs boson couplings
with SM gauge bosons, as discussed in Abdallah:2014fra  and seen in Fig. 1.

Figure 14:  Signal strength of the lightest and next-to-lightest Higgs bosons
in the BLSSM (in blue and red, respectively) in the $ZZ$ channel. The signal
strength of the lightest and next-to-lightest Higgs bosons in the MSSM are
given in violet and green points, respectively. The $1$ and $2\sigma$
confidence intervals are extracted from data collected during Run 1 with the
observed exclusion limit as given in Chatrchyan:2013mxa  is also included.

In Fig. 14, we show the signal strength of $h$ and $h^{\prime}$ decays to $ZZ$
for $m_{h}\approx 125$ GeV and $m_{h^{\prime}}$ around $140$ GeV along with 1
and $2\sigma$ confidence bands extracted from data collected during Run 1 with
the observed exclusion limit of Aad:2014fia . As the other two channels
previously discussed, the results of the BLSSM for both $h$ and $h^{\prime}$
match the observed data rather closely. We refrain from presenting here the
MSSM results for $h$ as in the decoupling limit they essentially coincide with
the SM ones (whereas those for the MSSM $H$ boson are outside the frame).

Figure 15: Number of signal events for $h$ and $h^{\prime}\to ZZ(\to 4\ell)$
decays ($\ell=e,\mu$) (red) induced by ggF and VBF versus the $4\ell$
invariant mass at $\sqrt{s}=13$ TeV after 100 fb-1 of luminosity alongside the
two dominant backgrounds (blue and black).

The results of our simulation for Run 2 are based on $ZZ\to 4\ell$ decays,
wherein $\ell=e,\mu$. In Fig. 15, we show the number of events for the $h$ and
$h^{\prime}$ bosons in the BLSSM plotted against the four-lepton invariant
mass. As can be seen from this plot, a promising signature of $h^{\prime}\to
ZZ\to 4\ell$ around $140\ {\rm GeV}$ emerges alongside the SM-like one at
$\approx 125$ GeV. The main contributions from SM backgrounds come from
$Z\gamma^{\ast}$ and $ZZ^{(*)}$. Significances at 100 fb-1 are already enough
to claim evidence in both Higgs channels.

Reconstruction of the $h$ and $h^{\prime}$ decays can only be performed for
one on-shell ($Z$) and one off-shell ($Z^{*}$) gauge boson, as
$M_{Z}<m_{h,h^{\prime}}<2M_{Z}$ for both Higgs states. We notice that the
combination of the two highest $p_{T}$ leptons is the most likely one to
emerge from the on-shell $Z$ boson decay while the other two leptons most
often come from the off-shell $Z$ boson decay. Fig. 16 shows the
reconstruction of both the off-shell and on-shell $Z$ boson decays for both
$h$ and $h^{\prime}$, illustrating that the off-shell distribution can be used
to increase the purity of each signal from cross-contamination.

In the light of such $Z$ boson spectra, we required the following cuts.

  1. 1.

The pseudorapidity of both electrons and muons is $|\eta|\leq 2.5$.

  2. 2.

We require a $Z$ candidate formed with a pair of leptons of the same flavour
and opposite charge, with mass window $40\leq M_{Z}\leq 120$ GeV, the
remaining leptons constructing the second off-shell $Z$ boson if they satisfy
$12\leq M_{Z}\leq 120$ GeV.

  3. 3.

In reconstructing the on-shell $Z$ we require the highest transverse momentum
lepton pair to be $\geq 20$ GeV.

  4. 4.

To protect the signals against leptons originating from hadron decays in jet
fragmentation or from the decay of low-mass hadronic resonances, we require
$M_{\ell^{+}\ell^{-}}\geq 4$ GeV, where $M_{\ell^{+}\ell^{-}}$ is the
invariant mass of any lepton pair.

Figure 16: Number of signal events for $h$ (red) and $h^{\prime}\to ZZ(\to
4\ell)$ (black) decays ($\ell=e,\mu$) induced by ggF and VBF versus the
$2\ell$ invariant mass at $\sqrt{s}=13$ TeV after 100 fb-1 of luminosity.
Left(Right): for the off(on)-shell $Z$ case.

Such a selection is already effective at 100 fb-1 and, as usual, increasing
luminosity will render this signal more and more significant, as per trend
seen in Fig. 17.

Figure 17: Left: Significance of the $h^{\prime}\to ZZ(\to 4\ell)$ signal (for
$m_{h^{\prime}}=140$ GeV and $\ell=e,\mu$) versus the luminosity (black).
Right: Number of events for signal and background for variable luminosity
(red). Data are produced at $\sqrt{s}=13$ TeV and the points correspond to an
integrated luminosity of 100, 300, 1000 and 3000 fb-1. Notice that event rates
are computed after the cuts described in the text.

##  V Conclusion

We have analysed the discovery potential of a second neutral Higgs boson in
the BLSSM at the LHC. We have confirmed that a double Higgs peak structure can
be accessed in this framework, in the $\gamma\gamma$,
$Z(\to\ell^{+}\ell^{-})\gamma$ and $ZZ(\to 4\ell)$ decay channels with Higgs
boson masses at $m_{h}\sim 125$ GeV and $m_{h^{\prime}}=140$ GeV, wherein $h$
and $h^{\prime}$ are the lightest CP-even Higgs states of the MSSM-like and
genuine BLSSM spectra, respectively.

Furthermore, under the assumption that the aforementioned excesses are not
confirmed by Run 2 data, we have studied the possibilities at the CERN machine
of establishing signals of an heavier $h^{\prime}$ state of the BLSSM. We have
shown that a peculiar decay in the BLSSM is $h^{\prime}\to hh$ (i.e., into a
pair of SM-like Higgs bosons), which can in fact be dominant from its
threshold (at $m_{h^{\prime}}\approx 2m_{h}\approx 250$ GeV) onwards. We have
shown that the associate $\gamma\gamma b\bar{b}$ signature can be
spectacularly visible over a wide mass interval, from, say, 250 to 500 GeV.

Combining all these results, and noting that similar Higgs signals would not
be available in the MSSM, we conclude that their extraction, either around 140
GeV or anywhere beyond 250 GeV or so, would not only point to a non-minimal
SUSY scenario, hence beyond the MSSM, but also possibly pinpoint the BLSSM.

Acknowledgements

SK is partially supported by the STDF project 13858. SM is supported in part
through the NExT Institute. All authors acknowledge support from the grant
H2020-MSCA-RISE-2014 n. 645722 (NonMinimalHiggs). AH is partially supported by
the EENP2 FP7-PEOPLE-2012-IRSES grant.

## References

  * (1) W. Abdallah, S. Khalil and S. Moretti, Phys. Rev. D 91, 014001 (2015).
  * (2) S. Khalil and S. Moretti, arXiv:1510.05934 [hep-ex].
  * (3) A. Hammad, S. Khalil and S. Moretti, arXiv:1503.05408 [hep-ph].
  * (4) S. Chatrchyan et al. [CMS Collaboration], Phys. Rev. D 89 (2014) 9, 092007.
  * (5) CMS Collaboration, CMS-PAS-HIG-13-001 (2013) and CMS-PAS-HIG-13-016 (2013).
  * (6) S. Chatrchyan et al. [CMS Collaboration], Phys. Lett. B 726, 587 (2013).
  * (7) S. Khalil and A. Masiero, Phys. Lett. B 665, 374 (2008).
  * (8) F. Staub, Adv. High Energy Phys. 2015 (2015) 840780
  * (9) F. Staub, Comput. Phys. Commun. 185, 1773 (2014).
  * (10) W. Porod and F. Staub, Comput. Phys. Commun. 183, 2458 (2012).
  * (11) J. S. Lee, M. Carena, J. Ellis, A. Pilaftsis and C. E. M. Wagner, Comput. Phys. Commun. 184, 1220 (2013).
  * (12) J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer, JHEP 1106, 128 (2011).
  * (13) T. Sjostrand, S. Mrenna and P. Z. Skands, Comput. Phys. Commun. 178, 852 (2008).
  * (14) J. de Favereau et al. [DELPHES 3 Collaboration], JHEP 1402, 057 (2014).
  * (15) E. Conte, B. Fuks and G. Serret, Comput. Phys. Commun. 184, 222 (2013).
  * (16) A. Belyaev, S. Khalil, S. Moretti and M. C. Thomas, JHEP 1405, 076 (2014).
  * (17) M. Hemeda, S. Khalil and S. Moretti, Phys. Rev. D 89, 011701 (2014).
  * (18) A. Djouadi, Phys. Rept. 459, 1 (2008).
  * (19) G. Aad et al. [ATLAS Collaboration], Phys. Rev. D 90, 112015 (2014).
  * (20) A. Djouadi, V. Driesen, W. Hollik and A. Kraft, Eur. Phys. J. C 1, 163 (1998).
  * (21) G. Aad et al. [ATLAS Collaboration], Phys. Lett. B 732, 8 (2014).
  * (22) ATLAS collaboration, ATLAS-CONF-2014-005.
  * (23) CMS Collaboration, CMS-PAS-HIG-14-013.
  * (24) G. Aad et al. [ATLAS Collaboration], Phys. Rev. Lett. 114, 081802 (2015).
  * (25) G. Aad et al. [ATLAS Collaboration], Eur. Phys. J. C 75, 412 (2015).
  * (26) CMS Collaboration, CMS-PAS-HIG-13-032.