aainstitutetext: School of Physical Sciences, Indian Association for the
Cultivation of Science,
2A & B Raja S.C. Mullick Road, Jadavpur, Kolkata 700 032,
Indiabbinstitutetext: Harish-Chandra Research Institute, A CI of Homi Bhabha
National Institute, Chhatnag Road, Jhunsi, Prayagraj (Allahabad) 211019,
Indiaccinstitutetext: Department of Theoretical Physics, Tata Institute of
Fundamental Research, 1 Homi Bhabha Road, Colaba, Mumbai 400005,
Indiaddinstitutetext: School of Physics and Institute for Collider Particle
Physics, University of the Witwatersrand, Johannesburg, Wits 2050, South
Africa

# Associated production of heavy Higgs bosons with a $b\bar{b}$ pair in the
Nonholomorphic MSSM and LHC searches

Utpal Chattopadhyay b AseshKrishna Datta c Samadrita Mukherjee d Abhaya Kumar
Swain<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In the NonHolomorphic Supersymmetric Standard Model (NHSSM), the Yukawa
couplings of the bottom quark ($y_{b}$) and the tau lepton ($y_{\tau}$) might
receive substantial supersymmetric (SUSY) radiative corrections which have
prominent dependencies on the NHSSM-specific trilinear soft parameters,
$A_{b}^{\prime}$ and $A_{\tau}^{\prime}$, respectively, in addition to their
well-known dependence on $\tan\beta$ as is already present in the Minimal SUSY
Standard Model (MSSM). We study to what extent these could affect the
production cross sections of the heavy Higgs bosons ($H$ and $A$) in
association with a pair of $b$-quarks and their decay branching ratios, in
particular, to a $\tau\bar{\tau}$ pair and compare them with those obtained in
the MSSM. Requiring compliance with the recently observed upper bounds on the
product of their total cross section and the branching ratio to
$\tau\bar{\tau}$ at the 13 TeV run of the Large Hadron Collider (LHC), with
data worth 139 fb-1, results in an altered exclusion region in the customary
$m_{A}-\tan\beta$ plane in the framework of the NHSSM when compared to what is
derived by the LHC experiments within an MSSM setup. Such alterations are
estimated to be pronounced only for large $\tan\beta$ ($\geq 40$) and for
large, negative $A_{b}^{\prime}$ when one finds a reinforced exclusion of the
$m_{A}-\tan\beta$ plane with an excluded $m_{A}$ value larger by $\approx 200$
GeV, compared to the MSSM case, at $\tan\beta=60$. On the other hand, the
maximum relaxation in $m_{A}$, for a similarly large but positive
$A_{b}^{\prime}$, barely exceeds $\approx 100$ GeV as a result of
complementary variations in production cross sections and decay branching
fractions of the heavy, neutral Higgs bosons.

###### Keywords:

Beyond Standard Model, Supersymmetry Phenomenology, Collider Physics

##  1 Introduction

Ever since its discovery at the Large Hadron Collider (LHC) the properties of
‘the’ Higgs boson are being probed in minute details thanks to the increasing
volume of the accumulated data. While its mass, which is not predicted by the
theory, is being estimated with an ever-increasing precision, various
couplings of the Higgs boson are being measured with even greater accuracy. It
is now known that whereas its couplings to the gauge bosons agree with their
predictions from the Standard Model (SM) at the $\sim 5$% level, the same with
the fermions from the third generation, to date, could attain compatibility
with the SM only at the level of 15-20% Aad:2019mbh .

The latter set of couplings are the Yukawa couplings ($y_{b}$, $y_{t}$ and
$y_{\tau}$) that are proportional to the respective fermion (bottom, top and
tau) masses. These are rather special in the sense that, when compared to the
minuscule ones for the light fermions from the first two generations, these
are much larger. In addition, the magnitudes of such couplings could also be
sensitive to scenarios beyond the SM, the implications of which could well
transcend the exclusive domain of particle physics by possibly influencing the
physics of the early Universe. Direct experimental probes to such couplings
thus have to be of paramount importance.

It is appreciated that any such fundamental coupling of a given excitation,
when could appear both in its production and decay, should be determined by
studying them in tandem to reduce the involved uncertainties in its
estimation. In the current context of the SM(-like) Higgs boson
($h_{\mathrm{SM}}$), such a combined approach is inherently unavoidable given
its production cross section at the LHC in any relevant mode and its decays
(branching fractions) could get to be largely interlaced and hence the
couplings are not independently extricable Wiesemann:2014ioa . As for the
Yukawa couplings, each of $y_{t}$ and $y_{b}$ might play a role in both
production and decay of the Higgs boson. However, as pointed out in reference
Wiesemann:2014ioa , the possibility of an efficient extraction of $y_{b}$ in
the decay $h_{\mathrm{SM}}\to b\bar{b}$ is marred by several issues even
though this has the largest branching fraction. A priori, a better sensitivity
to $y_{b}$ has been expected in the associated $h_{\mathrm{SM}}b\bar{b}$
production mode notwithstanding the large irreducible background that it
attracts. This is since such a mode is much less encumbered with theoretical
assumptions when compared to working with the decay branching fraction
BR[$h_{\mathrm{SM}}\to b\bar{b}$] which depends on to what extent
$h_{\mathrm{SM}}$ decays to other modes and can be highly model-dependent.
However, while the LHC experiments are yet to take a dedicated look into the
$h_{\mathrm{SM}}b\bar{b}$ production process (which has a cross section
similar to that of $h_{\mathrm{SM}}t\bar{t}$ production that has already been
studied by the experiments), a recent study Pagani:2020rsg  has cast a shadow
on its prospects.

Under the circumstances, studying productions of various Higgs bosons,
possibly heavier than $h_{\mathrm{SM}}$, in association with a $b\bar{b}$-pair
in scenarios with extended Higgs sector assumes a heightened significance.
Extensions of the SM endowed with two Higgs doublets, of which the Minimal
Supersymmetric extension of the SM (MSSM) is a popular example, are in
immediate reference here. Such production modes could be the direct probes to
such scenarios. This is since in these scenarios $y_{b}$ could get
significantly enhanced for large values of $\tan\beta\,(={v_{u}\over v_{d}})$,
the ratio of the vacuum expectation values (vev’s), $v_{u}$ and $v_{d}$, of
the neutral components of the two $SU(2)$ Higgs doublets with opposite
hypercharges, $h_{u}$ and $h_{d}$. The former gives tree-level masses to the
up-type quarks while the latter is responsible for the same for the down-type
quarks and the charged leptons of the scenario. In addition, the interaction
vertices of Higgs bosons and a $b\bar{b}$ pair have an additional
$1/\cos\beta$ dependence which lead to their enhanced interaction strengths
for larger values of $\tan\beta$. Even then, it has been observed that unless
$y_{b}$ gets significantly enhanced in such new physics scenarios, a
production process like $h_{\mathrm{SM}}b\bar{b}$ is unlikely to be useful for
the purpose. It has been recently pointed out Pagani:2020rsg  that such an
enhancement in $y_{b}$ is strongly disfavored by the current measurements of
the $h_{\mathrm{SM}}\to b\bar{b}$ decay. This is consistent with the
implications of the recent bounds on the $m_{A}$-$\tan\beta$ plane Aad:2020zxo
; Aaboud:2017sjh ; Sirunyan:2018zut , $m_{A}$ being the mass of the $CP$-odd
Higgs boson of the MSSM scenario, which indicates that the
$h_{\mathrm{SM}}b\bar{b}$ coupling may not be reinforced enough even when
$y_{b}$ is significantly enhanced because of a suppressed mixing between the
two $CP$-even Higgs states of the MSSM.

In this backdrop, a renewed impetus for carrying out such studies can now be
drawn from the recent discourses of an otherwise MSSM scenario but for the
presence of additional ‘nonholomorphic’ soft terms in its Lagrangian which
goes by the name of Nonholomorphic Supersymmetric Standard Model
(NHSSM)Girardello:1981wz ; Martin:1999hc ; Ross:2016pml ; Ross:2017kjc ;
Chattopadhyay:2017qvh ; Chakraborty:2019wav ; Haber:2007dj ; Bagger:1993ji ;
Ellwanger:1983mg ; Jack:1999ud . It is further known that the presence of
generic nonholomorphic terms like $A_{f}^{\prime}\phi^{*}\phi$ could have
important implications Chattopadhyay:2016ivr . In particular, the
nonholomorphic trilinear soft term from the bottom quark sector involving
$A_{b}^{\prime}$ and $A_{\tau}^{\prime}$ might play important roles in the
phenomenology of the NHSSM as these could radiatively alter the magnitudes of
$y_{b}$ and $y_{\tau}$, respectively Chattopadhyay:2018tqv , over and above
what is induced by $\tan\beta$ in the MSSM Hall:1993gn ; Hempfling:1993kv ;
Carena:1994bv ; Pierce:1996zz ; Logan:2000cz ; Antusch:2008tf  (see references
Carena:2002es  and Djouadi:2005gj  for an overview in the MSSM case). Even
then, the couplings $h_{\mathrm{SM}}b\bar{b}/\tau\bar{\tau}$ remain suppressed
because of the reason mentioned in the previous paragraph. In a complementary
way, the heavier Higgs bosons of the scenario are well-poised to receive the
maximal benefit of the suppressed Higgs mixing angle and hence their
productions in association with a $b\bar{b}$ pair might emerge to be more
potential probes to $y_{b}$. However, the LHC bound mentioned earlier that is
obtained assuming the MSSM framework is expected to be sensitive to the
parameter(s) (like $A_{b}^{\prime}$) of the NHSSM. Hence a recast of the bound
in the NHSSM scenario would be necessary before drawing a conclusion.

Radiatively corrected $y_{b}$ and $y_{\tau}$ in the NHSSM scenario are
expected to affect the constraints derived in the MSSM framework by the LHC
experiments on the heavy Higgs boson sector in the form of an exclusion in the
$m_{A}-\tan\beta$ plane. This is since these studies search for heavier Higgs
bosons in their productions in association with $b$-quark(s), i.e., $pp\to
b\bar{b}H/A$ Aad:2020zxo ; Aaboud:2017sjh ; Sirunyan:2018zut , $btH^{\pm}$
ATLAS:2018gfm ; ATLAS:2021upq , and their subsequent decays to $\tau$-leptons
(and to top quarks, in the case of $H^{\pm}$), i.e.,
$H,A\to\tau\bar{\tau},\,H^{\pm}\to\tau\nu_{\tau}$ ATLAS:2018gfm , $tb$
ATLAS:2021upq . However, an understanding of the extents to which $y_{b}$ and
$y_{\tau}$ could actually influence the proceedings requires a somewhat
detailed analysis which we take up duly in this work. As for the $H^{\pm}$
states, their searches at the LHC, as yet, could only set much weaker
exclusions on the $m_{A}-\tan\beta$ plane when compared to those from the
searches of the $H/A$, except for rather low values of $\tan\beta$ ($\lesssim
2$) which are anyway not so relevant for the phenomenology of the latter
states that we are interested in. Hence we do not discuss the related
phenomenology of the $H^{\pm}$ states separately in the present work.

Production processes of Higgs bosons in association with bottom squarks
($\tilde{b}$) also involve $y_{b}$ and hence are, a priori, in context.
However, we find that given the current ATLAS lower bounds on the mass of the
lighter sbottom ($\tilde{b}_{1}$) ATLAS:2021yij , the maximum
$\sigma_{(pp\to\tilde{b}_{1}{\tilde{b}}^{*}_{1}H)}$ obtainable over the NHSSM
parameter space for a relatively light ‘$H$’ state is smaller by quite a few
orders when compared to $\sigma_{(pp\to b\bar{b}H)}$ for the same $m_{H}$. For
our purposes, we add to the latter the contribution $\sigma_{(pp\to
b\bar{b}A)}$ where $m_{A}\simeq m_{H}$. In contrast, note that
$\sigma_{(pp\to\tilde{b}_{1}{\tilde{b}}^{*}_{1}A)}$ vanishes since the vertex
$A\tilde{b}_{1}\tilde{b}_{1}$ is absent Dedes:1998yt  and
$\sigma_{(pp\to\tilde{b}_{1}{\tilde{b}}^{*}_{2}A+\mathrm{c.c.})}$ is further
suppressed due to a more massive $\tilde{b}_{2}$ in the final state. Thus, the
combined rate for $H,A$ production in association with a pair of sbottoms is
way too small to be sensitive to the contemporary runs of the LHC. Hence we do
not consider those processes further. We are then set free to consider much
larger soft SUSY breaking masses for the sbottoms (as well as, for the staus),
$m_{\tilde{b},\tilde{\tau}}$. This, in turn, allows us to study a wider range
of values for the phenomenologically interesting NHSSM parameters like
$A_{b,\tau}^{\prime}\,(\lesssim m_{\tilde{b},\tilde{\tau}})$ without running
the risk of encountering a tachyonic spectrum or a charge and color breaking
(CCB) minimum in the scalar potential of the theory or jeopardizing the
stability of the electroweak vacuum Beuria:2017gtf . Henceforth, we refer to
these issues collectively as the CCB problem.

The paper is organized as follows. In section 2 we outline the theoretical
setup of the NHSSM scenario under study with an emphasis on how the same could
provide additional radiative corrections to running $y_{b,\tau}$ (in the
presence of a nonvanishing $A_{b,\tau}^{\prime}$) over and above the ones
found in the MSSM ($\tan\beta$-enhancement)) and how these new effects modify
the couplings of the neutral Higgs bosons to a $b$-quark or a $\tau$-lepton
pair. Section 3 is devoted to a detailed study of decays of the heavier
neutral Higgs bosons, ‘$H$’ and ‘$A$’, to a $b\bar{b}$ and a $\tau\bar{\tau}$
pair (the only accessible final states, following relevant experimental
studies) for varying $\tan\beta$ and $A_{b}^{\prime}$. This is followed by a
study of similar variations of the combined cross section for the associated
$b\bar{b}\Phi$ ($\Phi\equiv H,A$) productions at the 13 TeV LHC where we
discuss briefly some involved theoretical considerations that go into
obtaining a state-of-the-art, precise estimate of these cross sections. The
variation of the quantity $\sigma_{{}_{(pp\to
b\bar{b}H,A)}}\times\mathrm{BR}[H,A\to\tau\bar{\tau}]$ is then studied in
detail over the relevant parameter space of the NHSSM and then confronted by
the experimental bound on the same. The results are then recast on the
$m_{A}-\tan\beta$ plane with $A_{b,\tau}^{\prime}$ as parameters to find the
ways in which the exclusion contour obtained by the LHC experiments gets
modified in the NHSSM scenario. In section 4 we conclude.

##  2 The theoretical setup

In this section, we first briefly discuss the ingredients of the NHSSM
scenario which can be seen as an augmented version of the much popular MSSM
scenario. Given that the bottom and the tau Yukawa couplings are central to
the present study, we then outline how the radiative corrections to the
running masses (Yukawa couplings) $m_{b,\tau}$ ($y_{b,\tau}$) get modified in
the NHSSM scenario when compared to the corresponding MSSM contributions to
the same. Finally, we touch upon the structure of the Higgs sector which, by
itself, is no different from that in the MSSM scenario. Further, we point out
a few of its salient features that are in the direct context of the present
study.

###  2.1 The NHSSM scenario

The Lagrangian of the NHSSM scenario involving the soft SUSY breaking terms
includes, in addition to those that are already present in the MSSM,
nonholomorphic terms involving the scalars ($\phi^{2}\phi^{*}$) and the
fermions ($\psi\psi$). These are given by Martin:1999hc ;
Chattopadhyay:2016ivr

$\displaystyle-\mathcal{L^{\prime}}_{\text{soft}}^{\phi^{2}\phi^{*}}$
$\displaystyle=$ $\displaystyle\tilde{q}\cdot
h_{d}^{*}\,A_{u}^{\prime}\,\tilde{u}^{*}+\tilde{q}\cdot
h_{u}^{*}\,A_{d}^{\prime}\,\tilde{d}^{*}+\tilde{\ell}\cdot
h_{u}^{*}\,A_{\ell}^{\prime}\,\tilde{e}^{*}+\mathrm{h.c.}\;,$ (1)
$\displaystyle-\mathcal{L^{\prime}}_{\text{soft}}^{\psi\psi}$ $\displaystyle=$
$\displaystyle{\mathbf{\mu}^{\prime}}\tilde{h}_{u}\cdot\tilde{h}_{d}\;,$ (2)

where $\tilde{q}$, $\tilde{u}$ and $\tilde{d}$ represent the left-handed
(doublet) squark states, the right-handed (singlet) up- and down-type squark
states, respectively, while $\tilde{\ell}$ and $\tilde{e}$ stand for the left-
handed (doublet) and the right-handed (singlet) slepton states. $h_{u}$ and
$h_{d}$ are the Higgs fields as introduced in the previous section while
$\tilde{h}_{u}$ and $\tilde{h}_{d}$ are the corresponding higgsino states. The
trilinear parameters $A_{f=u,d,\ell}^{\prime}$ are defined as
$A^{\prime}_{f}=y_{f}\,{\cal A}^{\prime}_{f}$111The choice of the parameter
$A_{f}^{\prime}$ defined this way proves to be convenient as the numerical
packages that we employ later in our analysis adopt the same. This is true for
the trilinear parameters, $A_{f}$, of the MSSM as well, i.e.,
$A_{f}=y_{f}{\cal A}_{f}$. with ‘$f$’ denoting the type of the squark or the
slepton it is related to, $y_{f}$ and ${\cal A}^{\prime}_{f}$ being the Yukawa
coupling of the corresponding quark/lepton and its associated NH soft
trilinear parameter without scaling with the Yukawa coupling, respectively. It
should be noted that in contrast to the MSSM case, conjugate Higgs fields
$h_{d}^{*}$ and $h_{u}^{*}$ now appear in
$\mathcal{L^{\prime}}_{\text{soft}}^{\phi^{2}\phi^{*}}$ and to respect the
hypercharge assignments these now couple to up-type squarks and down-type
squarks/charged sleptons, respectively. $\mu^{\prime}$ represents the higgsino
mass parameter of the NHSSM which is a soft SUSY-breaking parameter of the
scenario, in contrast to its counterpart ‘$\mu$’ in the MSSM which appears in
the holomorphic, (SUSY-conserving) superpotential term
$\mu\tilde{h}_{u}\cdot\tilde{h}_{d}$.

The presence of the NH soft trilinear terms involving ${\cal A}^{\prime}_{f}$
modifies the off-diagonal terms of the tree-level sfermion mass-squared
matrices of the MSSM. Thus, in the NHSSM, these modified elements have the
generic form $-y_{f}\,v\left[{\cal A}_{f}-(\mu+{\cal
A}^{\prime}_{f})r_{\beta}\right]$ for the up- (with $r_{\beta}=\cot\beta$) and
the down-type squark/slepton (with $r_{\beta}=\tan\beta$) sectors
Chattopadhyay:2016ivr  with $v=\sqrt{v_{u}^{2}+v_{d}^{2}}=246$ GeV. Note that
unlike in the MSSM, contributions from the NH soft trilinear parameters ${\cal
A}^{\prime}_{f}$ arise through their products with either $\cot\beta$ (for the
up quark sector) or $\tan\beta$ (for the down quark and the charged lepton
sectors). Thus, for large $\tan\beta$, such contributions might be significant
in the bottom and the tau sectors. In turn, any phenomenological observable
that involves a chiral mixing of sfermions would be affected by such NH
parameters.

As for the neutralino and the chargino (the electroweakinos) sectors, those
are affected at the tree level in the presence of the NH $\mu^{\prime}$-term.
There, the modifications appear as shifts $\mu\to\mu+\mu^{\prime}$ in the
entries of the higgsino blocks of the respective mass matrices
Chattopadhyay:2016ivr .

###  2.2 Radiative corrections to $y_{b}$ and $y_{\tau}$

An important implication of the presence of NH soft trilinear terms involving
$A_{f}^{\prime}$ ($A_{b}^{\prime}$ and $A_{\tau}^{\prime}$) in the NHSSM
Lagrangian is that those could influence the radiative corrections to the
down-type Yukawa couplings, $y_{b}$ and $y_{\tau}$, thus modifying their tree-
level relations to $m_{b}$ and $m_{\tau}$, respectively. These are over and
above their dependencies on $\tan\beta$ that are already present in the MSSM.

At one-loop level, the principal MSSM contributions to $y_{b}$ arise from the
gluino-sbottom ($\tilde{g}-\tilde{b}_{i}$) loop and the chargino-stop
($\chi^{\pm}_{i}-\tilde{t}$) loop Hall:1993gn ; Hempfling:1993kv ;
Carena:1994bv ; Pierce:1996zz ; Logan:2000cz ; Antusch:2008tf .222The SUSY-QCD
and the SUSY-electroweak corrections can have enhanced contributions for
moderate to large ‘$\mu$’ and ${\cal A}_{t}$ when $\tan\beta$ is large. An
improved perturbative result is found by resumming these dominant
contributions which is also extended to include the effect from large ${\cal
A}_{b}$ that is absent at one-loop level Carena:1999py ; Guasch:2003cv . Given
that the parameters $\alpha_{{}_{S}}$ and $y_{t}$ that enter the corrections
have significant dependence on the renormalization scale, a further
improvement in the estimation with much reduced uncertainty is found by the
inclusion of two-loop SUSY-QCD corrections Noth:2008tw ; Noth:2010jy ;
Mihaila:2010mp ; Ghezzi:2017enb .  On the other hand, in the NHSSM where
$H_{u}$ is also associated with the down-type squarks (and sleptons) in the
trilinear soft interactions, one can have an extra $\tilde{g}-\tilde{b}$
contribution with a factor $y_{b}\,{\cal A}^{\prime}_{b}$ showing up at the
$\tilde{b}_{L}$-$H_{u}^{*}$-$\tilde{b}_{R}$ vertex. Besides, there is a
neutralino-sbottom ($\chi_{i}^{0}-\tilde{b}$) loop contribution with the same
factor $y_{b}\,{\cal A}^{\prime}_{b}$ appearing at the trilinear scalar vertex
and a factor of $y_{b}^{2}$ coming from the other two vertices having
$\tilde{b}_{L}\tilde{b}_{R}\chi^{0}$ interactions with their origins in the
superpotential Chattopadhyay:2018tqv . At one-loop level, the modified
relation between $m_{b}$ and $y_{b}$ in the NHSSM (which is identical in
appearance to that in the MSSM)333References Crivellin:2010er ;
Crivellin:2011jt ; Crivellin:2012zz  contain detailed analysis on the
effective Higgs boson vertices that involve nonholomorphic ${\cal
A}^{\prime}_{f}$ terms in SUSY-QCD and electroweak corrections to $m_{f}$
($y_{f}$). and hence the radiatively-corrected expression for $y_{b}$ are
given by

$m_{b}=\frac{y_{b}v_{d}}{\sqrt{2}}(1+\Delta_{b})\quad\Longrightarrow\quad
y_{b}=\frac{\sqrt{2}m_{b}}{v_{d}(1+\Delta_{b})}\,,$ (3)

where $\Delta_{b}$ is the SUSY threshold correction which now receives, in
addition to the MSSM contributions, the NHSSM ones and are to be evaluated at
the scale $M_{\text{SUSY}}$ given by the typical mass of the SUSY excitations
that run in the loops. For $\tan\beta>>1$ and $M_{\text{SUSY}}>>m_{Z}$,
$\Delta_{b}$ is given by444In the presence of an $\alpha_{s}$ driven
contribution and contributions that are enhanced for large values of $|{\cal
A}^{\prime}_{b}|$ and $|{\cal A}_{t}|$ (that we require to have the mass of
the SM-like Higgs boson in the right ballpark), we have dropped from the
expression of $\Delta_{b}$ the contributions proportional to the $U(1)$ and
$SU(2)$ gauge couplings, $g^{\prime}$ and $g$. However, all possible
contributions to such corrections are taken care of by the numerical package
we use for the purpose (see section 3.1).

$\displaystyle\Delta_{b}\,$
$\displaystyle\simeq\,\frac{y_{t_{\mathrm{SM}}}^{2}}{16\pi^{2}}\,\mu^{\prime\prime}{\cal
A}_{t}\,I(m_{\tilde{t}_{1}}^{2},m_{\tilde{t}_{2}}^{2},\mu^{\prime\prime
2})\tan\beta+\frac{2\alpha_{{}_{S}}}{3\pi}\,m_{\tilde{g}}\,(\mu+{\cal
A}^{\prime}_{b})\,I(m_{\tilde{b}_{1}}^{2},m_{\tilde{b}_{2}}^{2},m_{\tilde{g}}^{2})\tan\beta$
(4) $\displaystyle\hskip
152.0pt+\frac{y^{2}_{b_{\mathrm{SM}}}}{16\pi^{2}}\,\mu^{\prime\prime}\,(\mu+{\cal
A}^{\prime}_{b})\,I(m_{\tilde{b}_{1}}^{2},m_{\tilde{b}_{2}}^{2},\mu^{\prime\prime
2})\tan\beta\,,$

where $\mu^{\prime\prime}=\mu+\mu^{\prime}$ and the characteristic loop
(integral) function $I(a,b,c)$ involving loop-momentum ‘$u$’ and squared
masses $a,b,c$ of the three propagating states in the loop is given by
Hall:1993gn

$I(a,b,c)=\int_{0}^{\infty}\frac{udu}{(u+a)(u+b)(u+c)}=\frac{ab~{}\ln(a/b)+bc~{}\ln(b/c)+ca~{}\ln(c/a)}{(a-b)(b-c)(a-c)}\,,$
(5)

and $\alpha_{{}_{S}}$ is the strong coupling constant and ${\cal A}_{t}$ and
$m_{\tilde{g}}$ are the soft SUSY breaking holomorphic (of the MSSM) trilinear
coupling of the top squarks to the up-type Higgs field (without a scaling with
$y_{t}$) and the mass of the gluino, respectively. $m_{\tilde{t}_{1(2)}}$ and
$m_{\tilde{b}_{1(2)}}$ stand for the masses of the lighter (heavier) top and
bottom squark states, respectively. To make a generic numerical sense of this
involved function, it may be noted that the same can be parametrized as
$I(a,b,c)=K_{I}/a_{max}$ Carena:1994bv  where $a_{max}=\mathrm{max}(a,b,c)$ is
the maximum of the three squared masses appearing there and $K_{I}\simeq
0.5-0.9$ in the absence of a large hierarchy between these masses.555The
limiting values of $K_{I}=0.5\,(1)$ are obtained by carrying out the shown
integration with one (two) of the masses vanishing.

The corresponding corrections to $y_{\tau}$ ($m_{\tau}$) come from
$\chi^{\pm}_{i}$–$\tilde{\nu}_{\tau}$ and $\chi_{i}^{0}$–$\tilde{\tau}$ loops.
Similar to the case of $\Delta_{b}$, $\Delta_{\tau}$ is now driven by ${\cal
A}^{\prime}_{\tau}$ and $\tan\beta$ where, unlike in the former case, roles
played by $g^{\prime}$ and ‘$g$’ are important, in general. Thus, the
relations in equations 3 and 4 can be used for the modified expression for
$m_{\tau}$ ($y_{\tau}$) by dropping first the $\alpha_{{}_{S}}$-dependent
gluino contribution straightaway and replacing all bottom(top)-flavor related
parameters by their leptonic ($\tau$($\nu_{\tau}$)) counterparts before
including the $g^{\prime}$\- and $g$-dependent terms. As a result, the
corresponding first term which would represent the contribution from the
$\chi^{\pm}_{i}$–$\tilde{\nu}_{\tau}$ loop and would go as the (vanishing)
neutrino Yukawa coupling ($y_{t}\rightarrow y_{\nu_{\tau}}$) also drops out.
Tweaking appropriately the expressions for the $g^{\prime}$\- and
$g$-dependent MSSM contributions as presented in reference Girrbach:2009uy
for the NHSSM scenario, $\Delta_{\tau}$ takes the following form:

$\displaystyle\Delta_{\tau}$
$\displaystyle\simeq\frac{y^{2}_{\tau_{\mathrm{SM}}}}{16\pi^{2}}\,\mu^{\prime\prime}\,(\mu+{\cal
A}^{\prime}_{\tau})\,I(m_{\tilde{\tau}_{L}}^{2},m_{\tilde{\tau}_{R}}^{2},\mu^{\prime\prime
2})\tan\beta+\frac{g^{\prime 2}}{16\pi^{2}}\,\mu^{\prime\prime}\,(\mu+{\cal
A}^{\prime}_{\tau})\,I(m_{\tilde{\tau}_{L}}^{2},m_{\tilde{\tau}_{R}}^{2},\mu^{\prime\prime
2})\tan\beta$ $\displaystyle+\frac{g^{\prime
2}}{16\pi^{2}}\,M_{1}\,\left[\mu^{\prime\prime}\left\\{{1\over
2}I(M_{1}^{2},\mu^{\prime\prime
2},m_{\tilde{\tau}_{L}}^{2})-I(M_{1}^{2},\mu^{\prime\prime
2},m_{\tilde{\tau}_{R}}^{2})\right\\}+\mu\,I(M_{1}^{2},m_{\tilde{\tau}_{L}}^{2},m_{\tilde{\tau}_{R}}^{2})\right]\tan\beta$
$\displaystyle-\frac{g^{2}}{16\pi^{2}}\,M_{2}\,\mu^{\prime\prime}\,\left[{1\over
2}I(M_{2}^{2},\mu^{\prime\prime
2},m_{\tilde{\tau}_{L}}^{2})+I(M_{2}^{2},\mu^{\prime\prime
2},m_{\tilde{\nu}_{\tau}}^{2})\right]\tan\beta\,.$ (6)

The modified relation between $m_{\tau}$ and $y_{\tau}$ and hence the
expression for corrected $y_{\tau}$ are given by

$m_{\tau}=\frac{y_{\tau}v_{d}}{\sqrt{2}}(1+\Delta_{\tau})\quad\Longrightarrow\quad
y_{\tau}=\frac{\sqrt{2}m_{\tau}}{v_{d}(1+\Delta_{\tau})}\,.$ (7)

As can be seen from equations 3 – 7, the one-loop corrected $y_{b(\tau)}$ not
only has an explicit direct dependence on $\tan\beta$ as is the case in the
MSSM but also has an additional dependence on ${\cal A}^{\prime}_{b}\,({\cal
A}^{\prime}_{\tau})$ (via the quantity $\mu+{\cal A}^{\prime}_{b}$ ($\mu+{\cal
A}^{\prime}_{\tau}$)) in the NHSSM which is roughly on the same footing as
$\tan\beta$ when ‘$\mu$’ is not too large. Furthermore, it is clear from the
expressions in equations 3 and 7 that one expects enhancements (suppressions)
in $y_{b,\tau}$ (over their usual variations with $\tan\beta$) for
$\Delta_{b,\tau}<0\,(>0)$.

###  2.3 The Higgs sector

As mentioned earlier, the Higgs sector of the NHSSM is structurally similar to
that of the MSSM having two complex Higgs doublets of opposite hypercharges
thus leading to five physical Higgs states once the electroweak symmetry gets
broken. These are the two $CP$-even Higgs states, ‘$h$’ and ‘$H$’, of which
the lighter one ($h$) is SM-like ($h\equiv h_{\mathrm{SM}}$); a $CP$-odd one,
‘$A$’, and two charged Higgs states, $H^{\pm}$.

At the tree level, the masses of these physical states and their couplings
with the SM fermions and the gauge bosons depend only on two free parameters.
These are conventionally chosen to be $\tan\beta$ and the mass of the physical
$CP$-odd neutral Higgs state, $m_{A}$ Djouadi:2005gj . In terms of these two
parameters, the masses of the other four Higgs states, at the tree level, are
given by

$\displaystyle m_{h,H}^{2}$ $\displaystyle=$ $\displaystyle{1\over
2}\bigg{[}m_{A}^{2}+m_{Z}^{2}\mp\sqrt{(m_{A}^{2}+m_{Z}^{2})-4m_{A}^{2}m_{Z}^{2}\cos^{2}2\beta}\bigg{]},\quad
m_{H^{\pm}}^{2}=m_{A}^{2}+m_{W^{\pm}}^{2}\;,$ (8)

where $m_{Z}$ and $m_{W^{\pm}}$ are the masses of the SM ‘$Z$’ and $W^{\pm}$
bosons, respectively. It follows that for $m_{A}>>m_{Z}$, $m_{H}\approx
m_{H^{\pm}}\approx m_{A}$. When the Higgs sector conserves the $CP$ symmetry,
the physical $CP$-even neutral Higgs states ($h,H$) are obtained from their
weak counterparts via the rotation Gunion:1984yn ; Gunion:1989we

$\begin{pmatrix}H\\\ h\end{pmatrix}=\begin{pmatrix}\cos\alpha\;&\sin\alpha\\\
-\sin\alpha\;&\cos\alpha\end{pmatrix}\begin{pmatrix}H_{1}^{0}\\\
H_{2}^{0}\end{pmatrix}\,,$ (9)

where the Higgs mixing angle ‘$\alpha$’ controls the interactions of these
physical Higgs states with the fermions, the gauge bosons, the sfermions and
the electroweakinos and is given by $m_{A}$ and $\tan\beta$ as

$\alpha=\frac{1}{2}\tan^{-1}\left(\tan
2\beta\,\frac{m_{A}^{2}+m_{Z}^{2}}{m_{A}^{2}-m_{Z}^{2}}\right)\quad\mathrm{with}\quad-\frac{\pi}{2}\leq\alpha\leq
0\ .$ (10)

The couplings of these Higgs states with the bottom quark, which are in direct
reference in the present work, are given at the tree level by their MSSM
expressions Djouadi:2005gj

$\displaystyle G_{hbb}$ $\displaystyle=-i{m_{b}\over
v}{\sin\alpha\over\cos\beta}\;,\qquad G_{Hbb}=i{m_{b}\over
v}{\cos\alpha\over\cos\beta}\;,\qquad G_{Abb}={m_{b}\over
v}\tan\beta\,\gamma_{5}\;.$ (11)

In the units of the tree-level SM $b$-quark Yukawa coupling
$y_{b}={im_{b}\over v}$, the first two couplings shown in equation 11 reduce
to

$\displaystyle g_{hbb}$ $\displaystyle=$
$\displaystyle-{\sin\alpha\over\cos\beta}=\sin(\beta-\alpha)-\tan\beta\cos(\beta-\alpha)\;,$
$\displaystyle g_{Hbb}$ $\displaystyle=$
$\displaystyle{\cos\alpha\over\cos\beta}=\cos(\beta-\alpha)+\tan\beta\sin(\beta-\alpha)\,.$
(12)

Thus, in the so-called decoupling limit ($m_{A}>>m_{Z}$ and $\tan\beta>>1$),
when $\cos(\beta-\alpha)\to 0$, $\sin(\beta-\alpha)\to 1$ and the factor
$\gamma_{5}$ is ignored,666For a massive $b$-quark, the process of
dimensionally regularizing the calculation is somewhat involved given that the
Dirac matrix $\gamma_{5}$ is a four-dimensional construct. The resulting
effects of a finite $m_{b}$ are small $\sim{\cal O}\left({m_{b}^{2}\over
m_{A}^{2}}\right)$ Dawson:2005vi . one finds

$g_{hbb}\to 1,\qquad g_{Hbb}\to\tan\beta,\qquad g_{Abb}=\tan\beta\;.$ (13)

The bottom line is that the $Hbb$ coupling in the decoupling limit and the
$Abb$ coupling are both enhanced by a factor of $\tan\beta$. On the other
hand, the $hbb$ coupling gets $\tan\beta$-enhanced in the non-decoupling limit
for which $m_{A}\gtrsim m_{Z}$. It suffices to mention here that, in contrast,
the coupling of the Higgs bosons to the top quark(s) are generically
$\tan\beta$-suppressed (i.e., $\cot\beta$-enhanced) and the couplings of ‘$h$’
and ‘$H$’ to top quarks have just the opposite dependencies on $\tan\beta$
when compared to the same to bottom quarks.

At this point, it is important to note that these couplings receive crucial
higher-order (SUSY-QCD and SUSY-electroweak) corrections whereby the reduced
Yukawa couplings factors with resummed contributions (as discussed in section
2.2) included, for the neutral Higgs states, are given by Carena:2002es ;
Carena:1999py ; Guasch:2003cv

$\displaystyle g_{hbb}$ $\displaystyle\simeq$
$\displaystyle-{\sin\bar{\alpha}\over\cos\beta}\,\Big{[}1-{\Delta_{b}\over{1+\Delta_{b}}}\,(1+\cot\bar{\alpha}\cot\beta)\Big{]}\;,$
$\displaystyle g_{Hbb}$ $\displaystyle\simeq$
$\displaystyle{\cos\bar{\alpha}\over\cos\beta}\,\Big{[}1-{\Delta_{b}\over{1+\Delta_{b}}}\,(1-\tan\bar{\alpha}\cot\beta)\Big{]}\;,$
$\displaystyle g_{Abb}$ $\displaystyle\simeq$
$\displaystyle\tan\beta\,\Big{[}1-{\Delta_{b}\over{1+\Delta_{b}}}\,{1\over\sin^{2}\beta}\Big{]}\;,$
(14)

where $\bar{\alpha}$ stands for the radiatively corrected Higgs mixing angle
which now includes possible NHSSM effects in addition to the ones of the MSSM
origin and $\Delta_{b}$ is as defined in equation 4.

An analogous $\tan\beta$-dependence follows for $y_{\tau}$ but for the absence
of an $\alpha_{s}$-driven gluino contribution to $\Delta_{\tau}$ (as given in
equation 6; when compared to equation $\Delta_{b}$ of 4). Thus, $y_{\tau}$
would be intrinsically smaller and so would be its range of variation about
its nominal MSSM value when $A_{\tau}^{\prime}$ is varied, as opposed to the
range of variation of $y_{b}$ under a varying $A_{b}^{\prime}$. This implies
that it is $A_{b}^{\prime}$, rather than $A_{\tau}^{\prime}$, that would
broadly dictate the decay branching fractions of $H/A$ to $\tau\bar{\tau}$,
the final state which the LHC experiments exploit in their search for these
Higgs states in their associated productions with a $b\bar{b}$ pair
Aad:2020zxo ; Aaboud:2017sjh ; Sirunyan:2018zut .

Armed with the above knowledge of the $(H/A)\,b\bar{b}$, $(H/A)\tau\bar{\tau}$
couplings in the NHSSM, we now take up the study of the production of these
Higgs states in association with a $b\bar{b}$ pair at the LHC and their
subsequent decays to $b\bar{b}$ and $\tau\bar{\tau}$ pairs. The goal is to
find out how the variation of the quantity $\sigma_{(pp\to
b\bar{b}H/A)}\times\mathrm{BR}[H/A\to\tau\bar{\tau}]$ that offers the most
stringent LHC exclusion of the customary $m_{A}-\tan\beta$ plane differs in
the NHSSM from that in the MSSM in the presence of a non-vanishing
$A_{b}^{\prime}$ and hence the impact of the same in modifying the said
exclusion region. In the subsequent sections of this work where we present and
discuss our results, unless otherwise mentioned, we would refer to the
trilinear parameters $A_{b}^{\prime}$ and $A_{\tau}^{\prime}$ as defined in
section 2.

##  3 Results

In this section, we present our results by first studying the dependencies of
the Yukawa couplings $y_{b}$ and $y_{\tau}$ on the NHSSM trilinear parameters
$A_{b}^{\prime}$ and $A_{\tau}^{\prime}$, respectively, and $\tan\beta$. While
for the productions of ‘$H$’ and ‘$A$’ in association with a $b\bar{b}$ pair,
variations of both $y_{b}$ and $\tan\beta$ play major roles, their decays, are
expected to be further governed by $y_{\tau}$. A study on the decays of the
heavy Higgs bosons of the NHSSM follows. We then touch upon some important
theoretical issues in the calculation of the cross section for the process
$pp\to b\bar{b}H/A$ at the LHC energies before describing our approach of
calculating the same. All these are aimed to find out how the nontrivial
dependencies of various (Yukawa) couplings of such Higgs bosons affect the
sensitivity of the LHC to these states when produced in association with a
pair of bottom quarks. This then prompts us to examine to what extent the
relevant constraints obtained from the latest LHC searches in the said mode,
for an MSSM scenario, could get relaxed or further strengthened in the NHSSM
scenario of our interest. For the numerical analysis, we employ the publicly
available implementation of the NHSSM model777Some useful patches are kindly
provided by the authors of SARAH. in SARAH-v4.14.4 Staub:2013tta ;
Staub:2015kfa -generated SPheno Porod:2011nf .888Note that SARAH-generated
SPheno follows the convention $A_{f,f^{\prime}}\equiv y_{f}{\cal
A}_{f,f^{\prime}}$ which is adopted by us and is discussed earlier in section
2.1. Unless otherwise mentioned, all subsequent figures make use of the fixed
values and the ranges of variation (provided at the scale
$M_{\text{SUSY}}\simeq 3$ TeV) of various SUSY input parameters of Table 1.
These conform to $122\,\mathrm{GeV}<m_{h_{\mathrm{SM}}}<128\,\mathrm{GeV}$ and
pass the constraints from HiggsBounds-v5.9.1 Bechtle:2020pkv .

To focus on the NHSSM-specific sector of the parameter space, in particular,
the part that involves the trilinear coupling parameters $A_{b}^{\prime}$ and
$A_{\tau}^{\prime}$ and which could affect the phenomenology of the bottom
quark and tau lepton, we fix all other SUSY parameters to their suitable
values such that various important bounds from the LHC experiments are
complied with and the observed properties of the SM-like Higgs bosons could be
ensured. We treat $\tan\beta$ as a free parameter since its interplay with
$A_{b}^{\prime}$ and $A_{\tau}^{\prime}$ is central to our present study.

Note that we fix ‘$\mu$’ to not too large a value of 500 GeV so that the
scenario remains somewhat ‘natural’. Given that the contributions to
$\Delta_{b,\tau}$ depend significantly on the quantity $\mu+{\cal
A}^{\prime}_{b,\tau}$, we would thus look out for relatively large
$|A_{b,\tau}^{\prime}|$ to find $\Delta_{b,\tau}$ to be substantially large.
Furthermore, we choose a somewhat large value of $\mu^{\prime}$ (1.5 TeV) such
that the higgsino-like electroweakinos have masses around
$\sim\mu+\mu^{\prime}=2$ TeV. Together with $M_{1}$ and $M_{2}$ set to 2 TeV,
these would ensure that the electroweakinos are heavy enough not only to pass
the current experimental bounds on them but also to render themselves
inaccessible for the ‘$H$’ and the ‘$A$’ states decaying to a pair of them.

Parameters | MSSM | NHSSM
---|---|---
$\mu$ (GeV) | 500
$M_{1},M_{2},M_{3}$ (TeV) | 2, 2, 3
$\tan\beta$ | [5 : 60]
$m_{A}$ (TeV) | [0.35 : 3]
$m_{\tilde{q}},m_{\tilde{\ell}}$ (TeV) | 3
$A_{t,b,\tau}$ (TeV) | -2.7, 0, 0
$\mu^{\prime}$ (TeV) | – | 1.5
$A_{t}^{\prime}$ (TeV) | – | 0
$A_{b}^{\prime},A_{\tau}^{\prime}$ (TeV) | – | [-3 : 3]
Table 1: Fixed values and ranges adopted for various input SUSY parameters
provided to SPheno at the scale $M_{\text{SUSY}}\simeq 3$ TeV.

###  3.1 Variations of $bb\Phi$ and $\tau\tau\Phi$ interaction strengths

The characteristic dependencies of $y_{b}$ and $y_{\tau}$ and the same for the
resulting strengths of the $\Phi bb$ and $\Phi\tau\tau$ interactions,
respectively, where $\Phi\in\\{H,A\\}$, are discussed in detail in the
previous section. In particular, note that the corrections to $y_{b,\tau}$ and
the very form of the interactions of these heavy Higgs bosons to a bottom
quark or a tau lepton pair in the decoupling regime each lends a separate
$\tan\beta$-dependence to these couplings (see equation 14). In addition, as
discussed in section 2, in the NHSSM, corrections to $m_{b,\tau}$
($y_{b,\tau}$) depend on $A_{b,\tau}^{\prime}$. In this subsection we study
these dependencies on $\tan\beta$ and
$A_{b,\tau}^{\prime}\,(=y_{b,\tau}\,{\cal A}^{\prime}_{b,\tau})$.

Figure 1: Variations of $y_{f=b,\tau}$ as functions of (a)
$A^{\prime}_{f=b,\tau}$, for two fixed values of $\tan\beta$ (10 and 40) and
of (b) $\tan\beta$, for three fixed values of $A^{\prime}_{f=b,\tau}$ ($-3$
TeV, 0 and 3 TeV). Variation of $y_{b}$ ($y_{\tau}$) in the
$A_{b}^{\prime}(A_{\tau}^{\prime})-\tan\beta$ plane is shown in plot (c)
((d)).

First, we discuss the dependence of $y_{b,\tau}$ ($m_{b,\tau}$) on
$A_{b,\tau}^{\prime}$ and $\tan\beta$ arising from radiative corrections to
$m_{b,\tau}$. As has been noted at the end of section 2.2, a positive
(negative) $\Delta_{b,\tau}$ turns the radiatively corrected values of
$y_{b,\tau}$ smaller (larger) than their nominal values. In both cases,
effects are amplified for larger values of $\tan\beta$. In figure 1 we
illustrate the variations of $y_{b,\tau}$ as functions of
$A^{\prime}_{b,\tau}$ for two different values of $\tan\beta$ (=10 and 40).
The corresponding MSSM values of $y_{b,\tau}$, for any given value of
$\tan\beta$, are closely the ones found at the point of intersection of the
particular curve (line) and the line erected vertically at
$A^{\prime}_{b,\tau}=0$. On the other hand, figure 1 shows variations $y_{b}$
and $y_{\tau}$ as functions of $\tan\beta$ for three different values
$A^{\prime}_{b,\tau}$ ($=-3,0,3$ TeV). Again, lines with
$A^{\prime}_{b,\tau}=0$ (in black) would closely correspond to the variations
of $y_{b,\tau}$ with $\tan\beta$ in the MSSM. These two plots reveal the
individual magnitudes of $y_{b}$ and $y_{\tau}$ as they vary and the ranges of
their variations. On both counts, $y_{b}$ wins over $y_{\tau}$, hands down,
except when $A_{b}^{\prime}$ has a large positive value (3 TeV, in the present
case). In the latter case, $y_{b}$ always remains smaller than $y_{\tau}$ for
all values of $\tan\beta$ (see the magenta curve) and their difference grows
with growing $\tan\beta$. A much stronger dependence of $y_{b}$ on $\tan\beta$
and $A_{b}^{\prime}$, as compared to that of $y_{\tau}$ on $\tan\beta$ and
$A_{\tau}^{\prime}$ can be straightaway traced back to the corresponding
dependencies of $\Delta_{b}$ and $\Delta_{\tau}$, respectively, on the
relevant pairs from these parameters. These two plots further reveal that
while $y_{b}$ is sensitive to both $A_{b}^{\prime}$ and $\tan\beta$ (the
sensitivity to the former being milder than that to the latter) and assumes
larger (smaller) values for large negative (positive) $A_{b}^{\prime}$ and
large $\tan\beta$, $y_{\tau}$ has an almost exclusive and direct dependence on
$\tan\beta$ only. Figure 1 (1) captures the same variations but in the plane
of $A_{b}^{\prime}(A_{\tau}^{\prime})-\tan\beta$.

The resulting variation of the squared value of the reduced coupling strength
for $Hbb\,(\tau\tau)$ interaction, i.e., $g^{2}_{Hbb\,(\tau\tau)}$ (see
expressions under equation 14), in the decoupling regime, as a function of
$\tan\beta$, is shown in figure 2 (2) for
$A_{b}^{\prime}(A_{\tau}^{\prime})=(-3,0,3)$ TeV. These quantities are
extracted from SPheno outputs using equation 14. Note that while we are
compelled to adopt a logarithmic vertical axis for figure 2 since
$g^{2}_{Hbb}$ takes off rather sharply with a growing $\tan\beta$ (thanks to a
large negative $A_{b}^{\prime}$ and a catalyzing role played by the
$\alpha_{{}_{S}}$-dependent contribution to $\Delta_{b}$), we stick to a
linear one in figure 2 to illustrate the corresponding variations in
$g^{2}_{H\tau\tau}$ which are much less sensitive to the signs and magnitudes
of $A_{\tau}^{\prime}$. Curves with $A_{b}^{\prime},A_{\tau}^{\prime}=0$
roughly illustrate such variations in the MSSM-like setups and are shown for
reference purposes. As has been noted in section 2.3, for all practical
purposes for this work, one can consider
$g^{2}_{Abb\,(\tau\tau)}=g^{2}_{Hbb\,(\tau\tau)}$. Hence we do not present a
separate set of plots for $g^{2}_{Abb\,(\tau\tau)}$.

Given the factor $g_{Hbb\,(\tau\tau)}$ scales the $Hbb\,(\tau\tau)$ coupling
in the SM to its value in the NHSSM, $g^{2}_{Hbb\,(\tau\tau)}$ gives us a
quick estimate of the factors by which the SM-like values of the associated
production cross section $\sigma_{(pp\to b\bar{b}H/A)}$ and the relevant decay
widths $\Gamma(H/A\to b\bar{b},\tau\bar{\tau})$ are to be multiplied to find
their corresponding values in the NHSSM. In this work we applied such a
scaling on the readily available values of ${\sigma_{(pp\to
b\bar{b}H/A)}}|_{\mathrm{SM}}$ as a function of $m_{H/A}$ (see section 3.3) to
find their corresponding values in the NHSSM. On the other hand, although such
scalings for $\Gamma(H/A\to b\bar{b},\tau\bar{\tau})$ get automatically taken
care of in the SARAH-SPheno framework that we use, a quantitative knowledge of
the magnitudes of $g^{2}_{H/A\,bb(\tau\tau)}$ and their nature of variation
over the $A_{b}^{\prime}(A_{\tau}^{\prime})-\tan\beta$ plane aid a better
understanding of the physics interplay in the study of the all important
quantity, i.e., $\sigma_{(pp\to
b\bar{b}H/A)}\times\mathrm{BR}[H/A\to\tau\bar{\tau}]$.

Figure 2:  Variations of the squared scaling factors $g^{2}_{Hbb}$ (left) and
$g^{2}_{H\tau\tau}$ (right) as functions of $\tan\beta$ and for three fixed
values ($-3$ TeV, 0 and 3 TeV) of $A_{b}^{\prime}$ (left) and
$A_{\tau}^{\prime}$ (right).

Figure 3:  Variations of $g^{2}_{Hbb}$ as functions of $m_{\tilde{b}}$ =
$m_{\tilde{b}_{L}}$ = $m_{\tilde{b}_{R}}$ and for three fixed values of
$A_{b}^{\prime}$ for each case of (a) $A_{b}^{\prime}<0$ and (b)
$A_{b}^{\prime}>0$.

As for how the variations of these couplings affect the LHC studies in
context, it may be recalled that the exclusion region in the $m_{A}-\tan\beta$
plane from the searches of neutral heavy Higgs bosons ATLAS:2018gfm ;
Sirunyan:2018zut ; Aad:2020zxo ; Aaboud:2017sjh  is shaped by the variations
of the production cross section $\sigma_{(pp\to b\bar{b}H/A)}$ and the
branching fraction BR[$H/A\to\tau\bar{\tau}$] as functions of $m_{A}$ and
$\tan\beta$. As far as NHSSM-specific input parameters are concerned, while
the former is controlled by $A_{b}^{\prime}$, the latter are governed
collectively by $A_{b}^{\prime}$ and $A_{\tau}^{\prime}$. Thus, depending upon
the combinations of values of $A_{b}^{\prime}$ and $A_{\tau}^{\prime}$, the
reported region of exclusion in the $m_{A}-\tan\beta$ plane in the MSSM
scenario would get altered (either squeezed or extended) when the NHSSM
scenario is in context. Quantitatively, however, one can perhaps foresee that
the largeness of $\Delta_{b}$ (thanks to the SUSY-QCD correction which is
missing in $\Delta_{\tau}$) would control the proceedings except when the
input parameters get to be tuned to the contrary. We discuss the decays of the
‘$H$’ and the ‘$A$’ states in section 3.2 and their productions in association
with a $b\bar{b}$ pair in section 3.3.

Before we leave this subsection, it would be in order to have a brief
discussion on how the other relevant input SUSY parameters, viz.,
$m_{\tilde{g}}$ and $m_{\tilde{b}_{1,2}}$ could affect the phenomenology.
These parameters enter the proceedings effectively through the dominant
$\alpha_{{}_{S}}$-dependent term in $\Delta_{b}$ (see equation 4). In figure 3
we present the variations of $g^{2}_{\Phi bb}$ as a function of
$m_{\tilde{b}}=m_{\tilde{b}_{1}}=m_{\tilde{b}_{2}}$ for some fixed negative
(figure 3) and positive (figure 3) values of $A_{b}^{\prime}$, three of them
in each case. For each chosen value of $A_{b}^{\prime}$ we ensure
$m_{\tilde{b}}\geq|A_{b}^{\prime}|$ such that appearance of CCB minima could
be broadly avoided. Bands for different values of $A_{b}^{\prime}$ appear due
to a varying $m_{\tilde{g}}$ ($3\,\mathrm{TeV}\leq m_{\tilde{g}}\,\leq
6\,\mathrm{TeV}$). The top (bottom) edge of each band corresponds to the
largest (smallest) value of $m_{\tilde{g}}$ used in the scan when
$A_{b}^{\prime}<0$. The situation is reversed for $A_{b}^{\prime}>0$. The
following observations can be made from the plots in figure 3.

  * •

$g^{2}_{\Phi bb}$ decreases (increases) with increasing $m_{\tilde{b}}$ for
$A_{b}^{\prime}<0\,(>0)$.

  * •

$g^{2}_{\Phi bb}$ increases (decreases) with increasing $m_{\tilde{g}}$ for
$A_{b}^{\prime}<0\,(>0)$.

  * •

The smallest value ($\approx 1300$) obtained for $g^{2}_{\Phi bb}$ with
$A_{b}^{\prime}<0$ is larger than the largest value ($\approx 1000$) found for
the same with $A_{b}^{\prime}>0$.

Thus, in our present scheme of analysis, larger values of $g^{2}_{\Phi bb}$,
for $A_{b}^{\prime}<0$, are obtained for a larger
$|A_{b}^{\prime}|=m_{\tilde{b}}$ and for a larger $m_{\tilde{g}}$. The
possible maximum for $g^{2}_{\Phi bb}$, however, tends to saturate for larger
values of $|A_{b}^{\prime}|=m_{\tilde{b}}$. On the other hand, for
$A_{b}^{\prime}>0$, where we ask how small $g^{2}_{\Phi bb}$ could get, we
find that its smallest possible value $\sim 100$ (see figure 3) is approached
again for $|A_{b}^{\prime}|=m_{\tilde{b}}$ and for the largest $m_{\tilde{g}}$
that we choose. On this occasion, this value of $g^{2}_{\Phi bb}$ is found to
be nearly independent of $A_{b}^{\prime}$. These findings would come in handy
when we discuss in section 3.4 how the LHC constraints on the
$m_{A}-\tan\beta$ plane obtained in an MSSM scenario could get altered in the
NHSSM scenario.

###  3.2 Decays of the heavier neutral Higgs bosons

In the decoupling regime, the couplings of the heavier neutral Higgs bosons to
gauge bosons ($V\equiv W^{\pm},Z$), i.e., $g_{HVV}\propto\cos(\beta-\alpha)\to
0$. This results in a generic insensitivity in their searches in the modes
$H\to VV$, $A\to hZ$ etc. Under the circumstances, $H,A\to t\bar{t}$ would
dominate when these Higgs states are sufficiently heavy and $\tan\beta$ is not
large. With increasing $\tan\beta$ and hence enhanced $y_{b,\tau}$, the decays
$H,A\to b\bar{b},\tau\bar{\tau}$ emerge as the only possibilities unless some
lighter SUSY states (in particular, the relatively light electroweakinos) are
present in the spectrum to which ‘$H$’ and ‘$A$’ could also decay to with
moderate branching fractions, a feature that is already present in the MSSM.
However, we aim for a simpler setup suited for our purpose where such decays
are kinematically forbidden. This is achieved for all $m_{H,A}\leq 3$ TeV by
setting $\mu^{\prime}=1.5$ TeV, $M_{1,2}=2$ TeV and $M_{3}=3$ TeV (see table
1).

In the NHSSM, decays of $H/A$ to final states involving $b,\tau$ get further
affected if $A_{b}^{\prime},A_{\tau}^{\prime}\neq 0$ since these control the
magnitudes of $y_{b}$ and $y_{\tau}$, respectively. As expected, this could
well modify the experimental sensitivities of searches of heavier Higgs bosons
in these final states. Consequently, the exclusion region in the
$m_{A}-\tan\beta$ plane as obtained for the MSSM scenario would get altered
for the NHSSM case. Here, it is important to point out that although the
decays of heavier Higgs bosons to bottom quarks dominate at larger $\tan\beta$
(and large negative $A_{b}^{\prime}$), final states involving $b$-jets are
prone to large jetty backgrounds from the SM. Hence searches for the heavy
Higgs states at the LHC are best carried out in the subdominant modes
involving the $\tau$’s. Note, however, that while $y_{b}$ (and not $y_{\tau}$)
controls $\sigma(pp\to b\bar{b}H/A)$, BR[$H/A\to\tau\bar{\tau}$] are
controlled by both $y_{b}$ and $y_{\tau}$. The highest sensitivity to the
$\tau\bar{\tau}$ final state is thus attained when the right balance in the
magnitudes of $y_{b}$ and $y_{\tau}$ is struck over the NHSSM parameter space.

Given that the magnitudes and the patterns of variations of branching
fractions to $b\bar{b}$ and $\tau\bar{\tau}$ of the ‘$H$’ and ‘$A$’ states are
expected (and found) to be very similar, a common representative illustration
of these variations for the two states would suffice. In addition, the
complementarity of BR[$H/A\to b\bar{b}$] and BR[$H/A\to\tau\bar{\tau}$] allows
us to discuss the proceedings in terms of any one of these quantities. The
obvious choice is then BR[$H/A\to\tau\bar{\tau}$] since the concerned
experimental analyses consider the decays $H/A\to\tau\bar{\tau}$. It may
further be noted that BR[$H/A\to\tau\bar{\tau}$] depends not only on
$\tan\beta$ and $A_{\tau}^{\prime}$ but also on $A_{b}^{\prime}$. This is
since $A_{b}^{\prime}$ influences the partial decay width for $H/A\to
b\bar{b}$ which, in turn, affects BR[$H/A\to\tau\bar{\tau}$] by altering the
total decay width of $H/A$ that enters the computation of the branching
fractions.

(a) (b)

Figure 4: Variations of BR[$H/A\to\tau\bar{\tau}$] in the
$A_{b}^{\prime}-A_{\tau}^{\prime}$ plane for (a) $\tan\beta=10$ and (b)
$\tan\beta=40$ and for $m_{H,A}=2.5$ TeV.

In figure 4 we present the variations of BR[$H/A\to\tau\bar{\tau}$] in the
$A_{b}^{\prime}-A_{\tau}^{\prime}$ plane for $\tan\beta=10$ (figure 4(a)) and
$\tan\beta=40$ (figure 4(b)) when $m_{H,A}$ are set at a value of 2.5 TeV. It
is seen that, irrespective of the value of $\tan\beta$, the dependence of
BR[$H/A\to\tau\bar{\tau}$] on $A_{\tau}^{\prime}$ is only rather mild. Thus,
it is $A_{b}^{\prime}$ which almost single-handedly dictates
BR[$H/A\to\tau\bar{\tau}$]. As discussed in section 2.3, this can be traced
back to the smallness of $\Delta_{\tau}$ when compared to $\Delta_{b}$ as the
former lacks the gluonic (strong) contribution. It may further be noted how
rapidly BR[$H/A\to\tau\bar{\tau}$] drops as $A_{b}^{\prime}$ decreases from
its large positive value. Hence, for any given $\tan\beta$, the only situation
when decays to $\tau\bar{\tau}$ could be comparable to (or may even dominate)
those to $b\bar{b}$ is when $A_{b}^{\prime}$ has a large enough positive value
thus rendering the couplings $H/Ab\bar{b}$ much smaller compared to those for
$H/A\tau\bar{\tau}$. As can be found from figure 4, under such a circumstance,
BR[$H,A\to\tau\bar{\tau}$] could reach $\sim 90$% when $\tan\beta$ is large
while BR[$H,A\to\tau\bar{\tau}$] dropping to a mere $\sim 10$%.

It is important here to note that, for any given integrated luminosity, for
the LHC experiments to be sensitive to searches for the $H/A$ states in the
mode $pp\to b\bar{b}H,A$ followed by $H,A\to\tau\bar{\tau}$, the product
$\sigma(pp\to b\bar{b}H/A)\,\times$ BR[$H/A\to\tau\bar{\tau}$] has to be
optimally large. With $A_{b}^{\prime}$ effectively controlling both these
quantities but that being in a contrasting manner, only an optimal value of
$A_{b}^{\prime}$ would maximize the said product, for any given set of other
input parameters.

###  3.3 Production of heavier neutral Higgs bosons ($H,A$) with a bottom
quark pair

Predicting the production cross section of a Higgs boson in association with a
pair of $b$-quarks at hadron colliders involves theoretical subtleties.
$m_{b}$ being small yet may not be ignored, there are two competing approaches
to calculate the hard, parton level cross section: the so-called 4-flavor
scheme (4FS) and the 5-flavor scheme (5FS). These entail different
perturbation expansions of the same physical observables, i.e.,
$pp\to\mathrm{Higgs}+jets$, which start contributing at different perturbation
orders in $\alpha_{s}$. Their close agreement is then expected only at
sufficiently high orders Dawson:2005vi ; Campbell:2004pu .

In the 4FS, the $b$-quark is taken to be massive. Hence the proton (or the
antiproton) parton distribution contains gluons and the four light
($\sim$massless) flavor quarks, $u,d,s,c$ and the corresponding antiquarks.
The relevant parton level processes at the LHC, at the lowest order (LO) in
perturbation theory, are then driven by gluon-fusion and quark-antiquark
annihilation, i.e.,

$gg,q\bar{q}\to b\bar{b}~{}H/A\;.$ (15)

A generic, representative set of Feynman diagrams for this (4FS LO) category
is shown in the upper panel of figure 5. At the LHC energies, the LO cross
section is largely dominated by the gluon fusion process. Note that the
$t$-channel gluon-fusion process in figure 5 leads to large logarithms
$\sim\ln(\mu_{F}^{2}/m_{b}^{2})$ on integrating over the phase space of the
final state $b$-quarks; $\mu_{F}\sim m_{{}_{\Phi\in H,A}}$ being the so-called
factorization scale. In the computation of an inclusive cross section of the
sort, such large logarithms (involving finite $m_{b}$) are effectively
resummed to all orders via the introduction of the $b$-quark parton
distribution function (PDF) Barnett:1987jw  in the proton with $m_{b}\to 0$.
This is the 5FS. A representative set of Feynman diagrams in this scheme, for
the inclusive process in reference, is shown in the lower panel of figure 5.

4FS LO

Figure 5: Representative 4FS LO Feynman diagrams contributing to exclusive
$b\bar{b}\Phi$ ($\Phi=h,H,A$) production in the 4FS (top panel) and the ones
for the corresponding inclusive process at LO (left), NLO (middle) and NNLO
(right) in the 5FS (bottom panel). See text for details.

5FS

The fully inclusive cross section in the 5FS would then be dominated by the LO
process $b\bar{b}\to\Phi$ Dicus:1988cx  in figure 5 and exceeds the
corresponding 4FS LO estimates. This can be understood by noting that the
dominant 4FS LO process (the $t$-channel one of figure 5) appears only at the
next-to-next-to-leading order (NNLO)-QCD in the 5FS, as shown in 5. It is
further observed that the 4FS (5FS) approach is better suited in the
asymptotic limit of small (large) Higgs boson masses when the mass of the
bottom quark cannot (can) be neglected as is the defining characteristic of
the scheme Harlander:2011aa .

Experiments studying $pp\to b\bar{b}\Phi$ necessarily have to tag $b$-jet(s).
An exclusive search like this requires a minimum $p_{T}^{(b)}$ which ensures
that the tagged $b$-jets do not have their origins in the $b$-quark PDF. Hence
one may rely exclusively on the 4FS for the LO rates. This further guarantees
that the Higgs boson to have been radiated off one such final state $b$-quark
produced in the hard scattering thus probing $y_{b}$ directly.

Eventually, bounds derived from the experiments on $\sigma_{(pp\to
bbH/A)}\times\mathrm{BR}[H/A\to\tau\bar{\tau}])$, and on the $m_{A}-\tan\beta$
plane thereon ATLAS:2017eiz ; CMS:2018rmh ; ATLAS:2020zms , refer to the
inclusive (total) cross sections. In either of these two schemes, more precise
estimates of these cross sections, which are stable against variations in the
unphysical factorization/renormalization ($\mu_{R}$) scales, can only be found
by including higher order corrections.

Given that the computation of the higher-order corrections is simpler in the
5FS (since $m_{b}\rightarrow 0$ and the LO process is a resonant one), cross
sections up to NNLO (in $\alpha_{{}_{S}}$) are known there for long
Dicus:1998hs ; Balazs:1998sb ; Harlander:2003ai . In contrast, rates only up
to next-to-leading order (NLO) are available in the 4FS Dittmaier:2003ej ;
Dawson:2003kb ; Wiesemann:2014ioa .999This calculation is recently improved to
contain NLO-electroweak corrections Pagani:2020rsg . At these available
orders, the cross sections obtained in the two schemes tend to show improved
agreements Campbell:2004pu  along with a much milder dependence on $\mu_{F}$
and $\mu_{R}$. In spite of this, 4FS NLO results are unable to capture
logarithmic terms beyond the first few while the 5FS NNLO ones fail to do so
for the so-called power-suppressed terms Wiesemann:2014ioa .101010Of recent,
next-to-next-to-next-to leading order ($\mathrm{N^{3}LO}$) Duhr:2019kwi  and
N3-leading log (N3LL) Ajjath:2019neu  estimates in the 5FS have become
available and are found to agree even more closely with the NLO-QCD one in the
4FS.

Capturing the best from the two different kinematic regimes requires a
systematic combination of these contributions Harlander:2011aa ;
Aivazis:1993pi ; Thorne:1997ga ; Cacciari:1998it ; Kramer:2000hn ; Tung:2001mv
; Thorne:2006qt . Here, we adopt an improved approach called FNOLL-B
Forte:2015hba ; Bonvini:2015pxa ; Bonvini:2016fgf ; Forte:2016sja  as also
done by a recent ATLAS analysis ATLAS:2020zms . This approach combines the
4FS-NLO accuracy with the accuracy of NNLO leading-log (NNLL) for the resummed
collinear logarithms of the 5FS.111111Recently, an even more accurate
prediction of the inclusive production cross section by matching 4FS-NLO
results with those from the 5FS N3LO ones has been made Duhr:2020kzd . Towards
this, we follow the recommendations of references
LHCHiggsCrossSectionWorkingGroup:2016ypw ; recommend . Hence we use the cross
sections for $pp\to b\bar{b}H$ quoted there as a function of $m_{H}$
Bonvini:2015pxa ; Bonvini:2016fgf  where the $Hbb$ coupling is taken to be SM-
like and the common factorization ($\mu_{F}$) and the renormalization
($\mu_{R}$) scales, denoted by $\mu_{{}_{FR}}=\mu_{F}=\mu_{R}$, is set to
$(m_{H}+2m_{b})/4$.121212Note that $\mu_{{}_{FR}}$ is much smaller than
$m_{H}$ itself and is known to lead to reasonable agreements between 4FS and
5FS estimates already at LO and NLO Harlander:2003ai ; Maltoni:2003pn ;
Boos:2003yi  with improved perturbative convergence Bonvini:2015pxa ;
Maltoni:2003pn ; Maltoni:2005wd ; Maltoni:2012pa ; Harlander:2015xur . These
rates are indicated to be cross-checked against the FNOLL-B matched results of
reference Forte:2016sja .

We obtain $\sigma(pp\to b\bar{b}H/A)|_{\mathrm{\mbox{\tiny NHSSM}}}$ by an
appropriate rescaling of $y_{b}$ Dittmaier:2006cz ; Dawson:2011pe ;
Dittmaier:2014sva , i.e., by multiplying $\sigma(pp\to
b\bar{b}h_{\mathrm{SM}})$ by the square of the factors $g_{b\bar{b}H/A}$ of
equation 14 that are now computed for the NHSSM. It should be noted that, over
the broad range of $m_{H/A}$ we consider, their mass-differences are found to
be smaller than the experimental resolution (except when $m_{A}\leq 400$ GeV).
Hence these states are treated as mass-degenerate ATLAS:2020zms . Thus, for
any given mass $m_{H/A}$, their contributions are just added up.

In figure 6 we present the variations of the combined cross sections for
$pp\to b\bar{b}H,A$ at the 13 TeV LHC in the $m_{A}$–$\tan\beta$ plane for (a)
$A_{b}^{\prime}=3$ TeV, (b) $A_{b}^{\prime}=0$ (the MSSM-like case) and (c)
$A_{b}^{\prime}=-3$ TeV. As can be seen from these plots, for large
$\tan\beta$ and larger $m_{A}$, the cross sections with $A_{b}^{\prime}=-3$
TeV could be up to a couple of orders larger in magnitude than that when
$A_{b}^{\prime}=3$ TeV. The plot in the middle ($A_{b}^{\prime}=0$) serves
closely as an MSSM-reference for the purpose in which, for any given values of
$m_{A}$ and $\tan\beta$, one finds intermediate values of cross sections when
compared to the former two cases. The patterns of variations of the cross
section are intimately related to the same for the respective variation of
$y_{b}$ (see figure 1), or for that matter, of $g^{2}_{b\bar{b}H/A}$ (see
figure 2) in the same parameter-plane and these shape the projected exclusion
regions of the parameter space of the NHSSM scenario as derived in section
3.4.

Figure 6: Variations of cross sections $\sigma_{(pp\to b\bar{b}H+b\bar{b}A)}$
(in pb) at the 13 TeV LHC in the $m_{A}-\tan\beta$ plane for (a)
$A_{b}^{\prime}=3$ TeV, (b) $A_{b}^{\prime}=0$ and (c) $A_{b}^{\prime}=-3$
TeV.

###  3.4 Constraining $m_{A}-\tan\beta$ plane: NHSSM vis-a-vis LHC exclusion
(in the MSSM)

As pointed out in section 3.1, the sensitivity of a search for a heavy Higgs
boson via its associated production with a $b$-quark pair, at the LHC, broadly
depends on the quantity $\sigma(pp\to b\bar{b}H/A)$ times the branching
fraction of $H/A$ into an optimally sensitive mode which happens to be the
$\tau\bar{\tau}$ mode. Recent LHC analyses ATLAS:2017eiz ; CMS:2018rmh ;
ATLAS:2020zms  have put constraints (at 95% confidence level (CL)) on the
maximum allowed values of the product $\sigma(pp\to
b\bar{b}\Phi)\times\mathrm{BR}[\Phi\to\tau\bar{\tau}]$ as a function of
$m_{\Phi}$, in a model-independent way. These are then translated131313This is
possible since the narrow-width approximation is broadly reliable as
$\Gamma_{\Phi}\over m_{\Phi}$ ($\Phi\ni H,A$) remains on the smaller side
(varying between 0.5% – 5%) over the region of the NHSSM parameter space
chosen for the present study that includes large values of both $\tan\beta$
and $|A_{b}^{\prime}|$. into a model-dependent (MSSM) exclusion contour (again
at 95% CL) in the plane of $m_{A}$ and $\tan\beta$, the sole unknown
parameters of the theoretical setup in reference on which the observable
quantity under consideration depends. Given that both $\sigma(pp\to
b\bar{b}H/A)$ and $\mathrm{BR}[H/A\to\tau\bar{\tau}]$ crucially depend on an
NHSSM-specific parameter like $A_{b}^{\prime}$, in addition to their usual
(MSSM-type) dependencies on $m_{A}$ and $\tan\beta$, the existing domain of
exclusion in the $m_{A}-\tan\beta$ plane for the MSSM-like limit (i.e.,
$A_{b}^{\prime}=0$) is bound to get altered in the NHSSM as $A_{b}^{\prime}$
varies. We work out such exclusion contours in the NHSSM scenario and contrast
those with the one from a recent ATLAS analysis ATLAS:2020zms  that provides
the strongest exclusion till date on the $m_{A}-\tan\beta$ plane.

In figure 7 we study the variations of $\sigma\times\mathrm{BR}$ as functions
of $m_{{}_{H/A}}$ for two discrete values of each of $\tan\beta$ (10 and 40)
and $A_{b}^{\prime}$ ($-3$ TeV and 3 TeV)141414As mentioned earlier, ‘$H$’ and
‘$A$’ states could practically be considered mass-degenerate over the
considered range of their masses. Hence, for the purpose in hand, for any
given mass $m_{A}=m_{H}$, their production cross sections are added up and
then multiplied by a branching fraction to $\tau\bar{\tau}$ which happens to
be essentially the same for the two Higgs states, in particular, in the
scenario under consideration where their decays to only the SM states are
kinematically allowed. and contrast those with its observed upper limits (at
95% confidence limit (CL)) as reported by the ATLAS collaboration (in figure
2(b) of reference ATLAS:2020zms ) and represented by the lower/left edge of
the grey region in figure 7. The sections of the blue and the red curves
falling within the grey region indicate the values of $m_{{}_{H/A}}$ that are
ruled out in the NHSSM for different combinations of $\tan\beta$ and
$A_{b}^{\prime}$. As is expected, the largest (smallest) of the values of
$m_{{}_{H,A}}$ is ruled out with the dashed-blue (solid-red) curves for which
$\tan\beta$ is large (small) and $A_{b}^{\prime}$ is negative (positive).

Figure 7: (a) Varying $\sigma\times\mathrm{BR}$ for different values of
$A_{b}^{\prime}$ and $\tan\beta$, as a function of $m_{\Phi=H,A}$, shown
against its observed upper limits at 95% CL as reported in figure 2(b) of
reference ATLAS:2020zms  and defined by the lower/left edge of the grey region
and (b) the results of recasting the observed ATLAS upper limits on
$\sigma\times\mathrm{BR}$ in figure 7 on the $m_{A}-\tan\beta$ plane for
various representative values of $A_{b}^{\prime}$. In figures 7 and 7, curves
representing the reinforced (relaxed) bounds are obtained using the
combination $A_{b,\tau}^{\prime}=-(+)3$ TeV.

The same analysis for the NHSSM is then ported to delineate the altered
excluded regions in the $m_{A}-\tan\beta$ plane vis-a-vis the one presented by
the mentioned ATLAS experiment ATLAS:2020zms  for a generic MSSM scenario
called the $M_{h}^{125}$ scenario Bagnaschi:2018ofa  whose defining features
are all preserved in the NHSSM scenario we adopt in this work. While
contrasting our findings with the experimental one, we limit ourselves to the
associated $b\bar{b}H,A$ productions whereas the experiment considers the
gluon fusion process as well for the production of these Higgs states. It may,
however, be noted that at larger values of $\tan\beta$, where we find the
NHSSM effects to be somewhat prominent, $pp\to b\bar{b}H,A$ play the dominant
role Aaboud:2017sjh . Hence, given the scope of the present work, our approach
suffices. The closely overlapping dashed contours in green and black in figure
7 represent our derivation of the exclusion region with $A_{b}^{\prime}=0$
(the MSSM limit for the concerned sector) and the same obtained for the
$M_{h}^{125}$-MSSM scenario by the ATLAS experiment (at 95% CL; in figure 2(c)
of reference ATLAS:2020zms ), respectively. In NHSSM scenarios with
$A_{b}^{\prime}\neq 0$, one finds a slightly relaxed exclusion (to the edge of
the purple region) for $A_{b}^{\prime}>0$ and for all values of $\tan\beta$
(when compared to the same obtained for the $M_{h}^{125}$-MSSM scenario) and a
reinforced one (to the edge of the pink region) for $A_{b}^{\prime}<0$ only
when $\tan\beta$ is on the larger side ($\geqslant 30$). Note that for the
curves showing maximal deviations in $\sigma\times\mathrm{BR}$ in figure 7 and
hence maximal alterations in the exclusion limits in figure 7, with respect to
the $A_{b}^{\prime}=0$ case, are obtained with combinations of
$A_{b,\tau}^{\prime}$ with maximal values we use for them (i.e.,
$|A_{b,\tau}^{\prime}|=3$ TeV) and setting $\mathrm{sign}(A_{\tau}^{\prime})$
to be the same as that on $A_{b}^{\prime}$ which broadly controls the nature
of such alterations. We have checked that all these findings also agree
closely with the exclusion obtained using HiggsBounds-v5.9.1 which now
incorporates the pertinent ATLAS analysis carried out in reference
ATLAS:2020zms .

The nature of the modified exclusion contours could be broadly reconciled by
referring to the same for the variations of BR($H/A\to\tau\bar{\tau}$) and the
cross section $\sigma_{(pp\to b\bar{b}H/A)}$ as illustrated in figures 4 and
6, respectively. Of particular interest is the asymmetric nature of the
possible relaxations and reinforcements in the exclusions in the NHSSM
scenario that are obtained using somewhat large (extremal) values of
$A_{b}^{\prime}$ and $\tan\beta$, when compared to the exclusions obtained for
the MSSM case. Note from figure 4 that BR($H/A\to\tau\bar{\tau}$) more or less
saturates as $A_{b}^{\prime}$ and $\tan\beta$ approach large values
irrespective of the sign on $A_{b}^{\prime}$. However, while for large,
positive $A_{b}^{\prime}$, BR($H/A\to\tau\bar{\tau}$) reaches a plateau at a
larger value compared to the MSSM case, the reverse is true for large,
negative $A_{b}^{\prime}$. On top of that, from figures 6 and 6 (the MSSM
case) one can find that the cross sections in these two scenarios for any
given $\tan\beta$ and $A_{b}^{\prime}$ are not too different. Hence the
product $\sigma\times\mathrm{BR}$ remains nearly as sensitive to the
experiment as for the MSSM case. This is behind a minor relaxation in $m_{A}$
($\lesssim 100$ GeV) in the exclusion even though the setup is tailored to
bring about the maximum possible exclusion. On the other hand, it can be
gleaned from figure 6 that for large, negative $A_{b}^{\prime}$, the
difference in the cross sections $\sigma_{(pp\to b\bar{b}H/A)}$ with that from
the MSSM scenario quickly grows for growing $m_{A}$ and $\tan\beta$. Though
BR($H/A\to\tau\bar{\tau}$) tends to saturate at a value lower than that in the
MSSM case, it is entirely because of the rapid growth in the cross section
with $\tan\beta$ that drives the sensitivity to the experiment and hence a
stronger exclusion of the $m_{A}-\tan\beta$ plane (by $\sim 200$ GeV in
$m_{A}$, for $\tan\beta=60$) for large, negative $A_{b}^{\prime}$. At a more
fundamental level, these features are perhaps best understood in terms of the
variations of $y_{b,\tau}(\tan\beta,A_{b}^{\prime})$ as shown in figure 1 and
then, of $g^{2}_{Hbb}(\tan\beta,A_{b}^{\prime})$ as shown in figures 2 and 2.

In this context, it is noteworthy that the extents of relaxation and
reinforcement of exclusions that are presented in figure 7 both refer to NHSSM
scenarios that are cut out for yielding such alterations to be nearly maximal.
The only way such modifications could get somewhat pronounced is by setting
$m_{\tilde{g}}$ at a value higher than what we have used, i.e., 3 TeV (in
fact, this is the smallest value from the range we have used to draw the plots
in figure 3. On these plots, at $|A_{b}^{\prime}|=m_{\tilde{b}}$, the widths
of the bands are the largest and we are referring to the lowest (highest)
points from these bands at $|A_{b}^{\prime}|=m_{\tilde{b}}$ for
$A_{b}^{\prime}<0\,(>0)$.

##  4 Conclusions

In this work, we have discussed in detail the nature of the dependencies of
each of the bottom and the tau Yukawa couplings simultaneously on the NHSSM-
specific trilinear soft parameters $A_{b}^{\prime}$ and $A_{\tau}^{\prime}$,
respectively, and on $\tan\beta$. The MSSM higgsino mass parameter is held at
not too large a value (500 GeV) thus aiding the scenario to remain somewhat
‘natural’. We have further undertaken a thorough investigation on how these
could affect, in the NHSSM scenario, the production cross sections of the
heavy neutral Higgs bosons ($H$ and $A$) in association with a $b$-quark pair
and their branching fractions to $\tau\bar{\tau}$, the two key ingredients
that dictate their reach at the current and future runs of the LHC. These are
demonstrated in terms of altered (relaxed or extended) exclusion regions in
the customary $m_{A}-\tan\beta$ plane and contrasted with the latest such
constraint from the LHC for an MSSM scenario, for some representative values
of $A_{b}^{\prime}$. The role of $A_{\tau}^{\prime}$ is found to be rather
limited while the same of the masses of the sbottom and the gluino can be
moderate.

A general finding is that such deviations are maximal when both $\tan\beta$
and the magnitude of $A_{b}^{\prime}$ are large. For an example, for
$\tan\beta=60$, while a maximum relaxation of $\approx 100$ GeV in $m_{A}$ is
observed in the NHSSM scenario for a large positive $A_{b}^{\prime}$ ($=3$
TeV) when in the MSSM it is excluded up to 2 TeV, a strengthening of the
exclusion by about 200 GeV with respect to the latter occurs for a large
negative $A_{b}^{\prime}$ ($=-3$ TeV). For smaller values of $\tan\beta$, the
role of $A_{b}^{\prime}$ is set to diminish since, in that regime,
$A_{b}^{\prime}$ affects $\sigma(pp\to b\bar{b}H,A)$ and
BR[$H,A\to\tau\bar{\tau}$] in a complementary way such that the all-important
quantity $\sigma\times\mathrm{BR}$ ceases to be sensitive to a variation in
$A_{b}^{\prime}$. Thus, for smaller values of $\tan\beta$, the exclusion
contour in the $m_{A}-\tan\beta$ plane in the NHSSM scenario does not deviate
much from the reported (MSSM) one. Overall, unless for a rather large value of
$\tan\beta$, the relevant NHSSM parameters tend to conspire such that the LHC
sensitivity of the scenario does not differ much from the MSSM case. As far as
the exclusion of the $m_{A}-\tan\beta$ parameter plane is concerned, this
pattern is expected to be broadly preserved at the upcoming runs of the LHC.

## Acknowledgements

AKS thanks the Department of Physics, SGTB Khalsa College, University of Delhi
and the SERB sponsored Multi Institutional Project titled “Probing New Physics
Interactions” (CRG/2018/004889) running there under which the major part of
the work has been carried out. The authors also like to thank Mark Goodsell
for many useful discussions on various aspects of the implementation of the
NHSSM scenario in SPheno via SARAH. They also thank Arnaud Ferrari from the
ATLAS Collaboration, Rikkert Frederix, Jean-Loic Kneur, Davide Napoletano,
Werner Porod, V. Ravindran and Michael Spira for very helpful exchanges.

## References

  * (1) G. Aad et al. [ATLAS], Phys. Rev. D 101 (2020) no.1, 012002 doi:10.1103/PhysRevD.101.012002 [arXiv:1909.02845 [hep-ex]].
  * (2) M. Wiesemann, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni and P. Torrielli, JHEP 02 (2015), 132 doi:10.1007/JHEP02(2015)132 [arXiv:1409.5301 [hep-ph]].
  * (3) D. Pagani, H. S. Shao and M. Zaro, JHEP 11 (2020), 036 doi:10.1007/JHEP11(2020)036 [arXiv:2005.10277 [hep-ph]].
  * (4) A. M. Sirunyan et al. [CMS], JHEP 09 (2018), 007 doi:10.1007/JHEP09(2018)007 [arXiv:1803.06553 [hep-ex]].
  * (5) G. Aad et al. [ATLAS], Phys. Rev. Lett. 125 (2020) no.5, 051801 doi:10.1103/PhysRevLett.125.051801 [arXiv:2002.12223 [hep-ex]].
  * (6) M. Aaboud et al. [ATLAS], JHEP 01 (2018), 055 doi:10.1007/JHEP01(2018)055 [arXiv:1709.07242 [hep-ex]].
  * (7) L. Girardello and M. T. Grisaru, Nucl. Phys. B 194 (1982), 65 doi:10.1016/0550-3213(82)90512-0
  * (8) S. P. Martin, Phys. Rev. D 61 (2000), 035004 doi:10.1103/PhysRevD.61.035004 [arXiv:hep-ph/9907550 [hep-ph]].
  * (9) G. G. Ross, K. Schmidt-Hoberg and F. Staub, Phys. Lett. B 759 (2016) 110 doi:10.1016/j.physletb.2016.05.053 [arXiv:1603.09347 [hep-ph]].
  * (10) G. G. Ross, K. Schmidt-Hoberg and F. Staub, JHEP 1703 (2017) 021 doi:10.1007/JHEP03(2017)021 [arXiv:1701.03480 [hep-ph]].
  * (11) U. Chattopadhyay, D. Das and S. Mukherjee, JHEP 01 (2018), 158 doi:10.1007/JHEP01(2018)158 [arXiv:1710.10120 [hep-ph]].
  * (12) S. Chakraborty and T. S. Roy, Phys. Rev. D 100, no.3, 035020 (2019) doi:10.1103/PhysRevD.100.035020 [arXiv:1904.10144 [hep-ph]].
  * (13) H. E. Haber and J. D. Mason, Phys. Rev. D 77 (2008), 115011 doi:10.1103/PhysRevD.77.115011 [arXiv:0711.2890 [hep-ph]].
  * (14) J. Bagger and E. Poppitz, Phys. Rev. Lett. 71 (1993), 2380-2382 doi:10.1103/PhysRevLett.71.2380 [arXiv:hep-ph/9307317 [hep-ph]].
  * (15) U. Ellwanger, Phys. Lett. B 133 (1983), 187-191 doi:10.1016/0370-2693(83)90557-9
  * (16) I. Jack and D. R. T. Jones, Phys. Lett. B 457 (1999), 101-108 doi:10.1016/S0370-2693(99)00530-4 [arXiv:hep-ph/9903365 [hep-ph]].
  * (17) U. Chattopadhyay and A. Dey, JHEP 10 (2016), 027 doi:10.1007/JHEP10(2016)027 [arXiv:1604.06367 [hep-ph]].
  * (18) U. Chattopadhyay, A. Datta, S. Mukherjee and A. K. Swain, JHEP 10 (2018), 202 doi:10.1007/JHEP10(2018)202 [arXiv:1809.05438 [hep-ph]].
  * (19) L. J. Hall, R. Rattazzi and U. Sarid, Phys. Rev. D 50 (1994), 7048-7065 doi:10.1103/PhysRevD.50.7048 [arXiv:hep-ph/9306309 [hep-ph]].
  * (20) R. Hempfling, Phys. Rev. D 49 (1994), 6168-6172 doi:10.1103/PhysRevD.49.6168
  * (21) M. Carena, M. Olechowski, S. Pokorski and C. E. M. Wagner, Nucl. Phys. B 426 (1994), 269-300 doi:10.1016/0550-3213(94)90313-1 [arXiv:hep-ph/9402253 [hep-ph]].
  * (22) D. M. Pierce, J. A. Bagger, K. T. Matchev and R. j. Zhang, Nucl. Phys. B 491 (1997), 3-67 doi:10.1016/S0550-3213(96)00683-9 [arXiv:hep-ph/9606211 [hep-ph]].
  * (23) H. E. Logan, Nucl. Phys. B Proc. Suppl. 101 (2001), 279-288 doi:10.1016/S0920-5632(01)01512-2 [arXiv:hep-ph/0102029 [hep-ph]].
  * (24) S. Antusch and M. Spinrath, Phys. Rev. D 78 (2008), 075020 doi:10.1103/PhysRevD.78.075020 [arXiv:0804.0717 [hep-ph]].
  * (25) M. Carena and H. E. Haber, Prog. Part. Nucl. Phys. 50, 63-152 (2003) doi:10.1016/S0146-6410(02)00177-1 [arXiv:hep-ph/0208209 [hep-ph]].
  * (26) A. Djouadi, Phys. Rept. 459 (2008), 1-241 doi:10.1016/j.physrep.2007.10.005 [arXiv:hep-ph/0503173 [hep-ph]].
  * (27) M. Aaboud et al. [ATLAS], JHEP 09 (2018), 139 doi:10.1007/JHEP09(2018)139 [arXiv:1807.07915 [hep-ex]].
  * (28) G. Aad et al. [ATLAS], JHEP 06 (2021), 145 doi:10.1007/JHEP06(2021)145 [arXiv:2102.10076 [hep-ex]].
  * (29) G. Aad et al. [ATLAS], JHEP 05 (2021), 093 doi:10.1007/JHEP05(2021)093 [arXiv:2101.12527 [hep-ex]].
  * (30) A. Dedes and S. Moretti, Phys. Rev. D 60 (1999), 015007 doi:10.1103/PhysRevD.60.015007 [arXiv:hep-ph/9812328 [hep-ph]].
  * (31) J. Beuria and A. Dey, JHEP 10 (2017), 154 doi:10.1007/JHEP10(2017)154 [arXiv:1708.08361 [hep-ph]].
  * (32) M. Carena, D. Garcia, U. Nierste and C. E. M. Wagner, Nucl. Phys. B 577, 88-120 (2000) doi:10.1016/S0550-3213(00)00146-2 [arXiv:hep-ph/9912516 [hep-ph]].
  * (33) J. Guasch, P. Hafliger and M. Spira, Phys. Rev. D 68, 115001 (2003) doi:10.1103/PhysRevD.68.115001 [arXiv:hep-ph/0305101 [hep-ph]].
  * (34) D. Noth and M. Spira, Phys. Rev. Lett. 101, 181801 (2008) doi:10.1103/PhysRevLett.101.181801 [arXiv:0808.0087 [hep-ph]].
  * (35) D. Noth and M. Spira, JHEP 06, 084 (2011) doi:10.1007/JHEP06(2011)084 [arXiv:1001.1935 [hep-ph]].
  * (36) M. Ghezzi, S. Glaus, D. Müller, T. Schmidt and M. Spira, Eur. Phys. J. C 81, no.3, 259 (2021) doi:10.1140/epjc/s10052-021-09035-6 [arXiv:1711.02555 [hep-ph]].
  * (37) L. Mihaila and C. Reisser, JHEP 08, 021 (2010) doi:10.1007/JHEP08(2010)021 [arXiv:1007.0693 [hep-ph]].
  * (38) A. Crivellin, Phys. Rev. D 83 (2011), 056001 doi:10.1103/PhysRevD.83.056001 [arXiv:1012.4840 [hep-ph]].
  * (39) A. Crivellin, L. Hofer and J. Rosiek, JHEP 07 (2011), 017 doi:10.1007/JHEP07(2011)017 [arXiv:1103.4272 [hep-ph]].
  * (40) A. Crivellin and C. Greub, Phys. Rev. D 87 (2013), 015013 [erratum: Phys. Rev. D 87 (2013), 079901] doi:10.1103/PhysRevD.87.015013 [arXiv:1210.7453 [hep-ph]].
  * (41) J. Girrbach, S. Mertens, U. Nierste and S. Wiesenfeldt, JHEP 05 (2010), 026 doi:10.1007/JHEP05(2010)026 [arXiv:0910.2663 [hep-ph]].
  * (42) J. F. Gunion and H. E. Haber, Nucl. Phys. B 272 (1986), 1 [erratum: Nucl. Phys. B 402 (1993), 567-569] doi:10.1016/0550-3213(86)90340-8
  * (43) J. F. Gunion, H. E. Haber, G. L. Kane and S. Dawson, Front. Phys. 80 (2000), 1-404 SCIPP-89/13.
  * (44) S. Dawson, C. B. Jackson, L. Reina and D. Wackeroth, Mod. Phys. Lett. A 21 (2006), 89-110 doi:10.1142/S0217732306019256 [arXiv:hep-ph/0508293 [hep-ph]].
  * (45) F. Staub, Comput. Phys. Commun. 185 (2014), 1773-1790 doi:10.1016/j.cpc.2014.02.018 [arXiv:1309.7223 [hep-ph]].
  * (46) F. Staub, Adv. High Energy Phys. 2015 (2015), 840780 doi:10.1155/2015/840780 [arXiv:1503.04200 [hep-ph]].
  * (47) W. Porod and F. Staub, Comput. Phys. Commun. 183 (2012), 2458-2469 doi:10.1016/j.cpc.2012.05.021 [arXiv:1104.1573 [hep-ph]].
  * (48) P. Bechtle, D. Dercks, S. Heinemeyer, T. Klingl, T. Stefaniak, G. Weiglein and J. Wittbrodt, Eur. Phys. J. C 80 (2020) no.12, 1211 doi:10.1140/epjc/s10052-020-08557-9 [arXiv:2006.06007 [hep-ph]].
  * (49) J. M. Campbell, S. Dawson, S. Dittmaier, C. Jackson, M. Kramer, F. Maltoni, L. Reina, M. Spira, D. Wackeroth and S. Willenbrock, [arXiv:hep-ph/0405302 [hep-ph]].
  * (50) R. M. Barnett, H. E. Haber and D. E. Soper, Nucl. Phys. B 306 (1988), 697-745 doi:10.1016/0550-3213(88)90440-3
  * (51) D. A. Dicus and S. Willenbrock, Phys. Rev. D 39 (1989), 751 doi:10.1103/PhysRevD.39.751
  * (52) R. Harlander, M. Kramer and M. Schumacher, [arXiv:1112.3478 [hep-ph]].
  * (53) M. Aaboud et al. [ATLAS], JHEP 01 (2018), 055 doi:10.1007/JHEP01(2018)055 [arXiv:1709.07242 [hep-ex]].
  * (54) A. M. Sirunyan et al. [CMS], JHEP 09 (2018), 007 doi:10.1007/JHEP09(2018)007 [arXiv:1803.06553 [hep-ex]].
  * (55) G. Aad et al. [ATLAS], Phys. Rev. Lett. 125 (2020) no.5, 051801 doi:10.1103/PhysRevLett.125.051801 [arXiv:2002.12223 [hep-ex]].
  * (56) D. Dicus, T. Stelzer, Z. Sullivan and S. Willenbrock, Phys. Rev. D 59 (1999), 094016 doi:10.1103/PhysRevD.59.094016 [arXiv:hep-ph/9811492 [hep-ph]].
  * (57) C. Balazs, H. J. He and C. P. Yuan, Phys. Rev. D 60 (1999), 114001 doi:10.1103/PhysRevD.60.114001 [arXiv:hep-ph/9812263 [hep-ph]].
  * (58) R. V. Harlander and W. B. Kilgore, Phys. Rev. D 68 (2003), 013001 doi:10.1103/PhysRevD.68.013001 [arXiv:hep-ph/0304035 [hep-ph]].
  * (59) S. Dittmaier, M. Krämer and M. Spira, Phys. Rev. D 70 (2004), 074010 doi:10.1103/PhysRevD.70.074010 [arXiv:hep-ph/0309204 [hep-ph]].
  * (60) S. Dawson, C. B. Jackson, L. Reina and D. Wackeroth, Phys. Rev. D 69 (2004), 074027 doi:10.1103/PhysRevD.69.074027 [arXiv:hep-ph/0311067 [hep-ph]].
  * (61) C. Duhr, F. Dulat and B. Mistlberger, Phys. Rev. Lett. 125 (2020) no.5, 051804 doi:10.1103/PhysRevLett.125.051804 [arXiv:1904.09990 [hep-ph]].
  * (62) A. H. Ajjath, A. Chakraborty, G. Das, P. Mukherjee and V. Ravindran, JHEP 11 (2019), 006 doi:10.1007/JHEP11(2019)006 [arXiv:1905.03771 [hep-ph]].
  * (63) M. A. G. Aivazis, J. C. Collins, F. I. Olness and W. K. Tung, Phys. Rev. D 50 (1994), 3102-3118 doi:10.1103/PhysRevD.50.3102 [arXiv:hep-ph/9312319 [hep-ph]].
  * (64) R. S. Thorne and R. G. Roberts, Phys. Rev. D 57 (1998), 6871-6898 doi:10.1103/PhysRevD.57.6871 [arXiv:hep-ph/9709442 [hep-ph]].
  * (65) M. Cacciari, M. Greco and P. Nason, JHEP 05 (1998), 007 doi:10.1088/1126-6708/1998/05/007 [arXiv:hep-ph/9803400 [hep-ph]].
  * (66) M. Krämer, F. I. Olness and D. E. Soper, Phys. Rev. D 62 (2000), 096007 doi:10.1103/PhysRevD.62.096007 [arXiv:hep-ph/0003035 [hep-ph]].
  * (67) W. K. Tung, S. Kretzer and C. Schmidt, J. Phys. G 28 (2002), 983-996 doi:10.1088/0954-3899/28/5/321 [arXiv:hep-ph/0110247 [hep-ph]].
  * (68) R. S. Thorne, Phys. Rev. D 73 (2006), 054019 doi:10.1103/PhysRevD.73.054019 [arXiv:hep-ph/0601245 [hep-ph]].
  * (69) S. Forte, D. Napoletano and M. Ubiali, Phys. Lett. B 751 (2015), 331-337 doi:10.1016/j.physletb.2015.10.051 [arXiv:1508.01529 [hep-ph]].
  * (70) M. Bonvini, A. S. Papanastasiou and F. J. Tackmann, JHEP 11 (2015), 196 doi:10.1007/JHEP11(2015)196 [arXiv:1508.03288 [hep-ph]].
  * (71) M. Bonvini, A. S. Papanastasiou and F. J. Tackmann, JHEP 10 (2016), 053 doi:10.1007/JHEP10(2016)053 [arXiv:1605.01733 [hep-ph]].
  * (72) S. Forte, D. Napoletano and M. Ubiali, Phys. Lett. B 763 (2016), 190-196 doi:10.1016/j.physletb.2016.10.040 [arXiv:1607.00389 [hep-ph]].
  * (73) C. Duhr, F. Dulat, V. Hirschi and B. Mistlberger, JHEP 08 (2020) no.08, 017 doi:10.1007/JHEP08(2020)017 [arXiv:2004.04752 [hep-ph]].
  * (74) D. de Florian et al. [LHC Higgs Cross Section Working Group], doi:10.23731/CYRM-2017-002 [arXiv:1610.07922 [hep-ph]].
  * (75) See ‘‘https://twiki.cern.ch/twiki/bin/view/LHCPhysics/LHCHWGBBH#NLO_NNLLpart_ybyt_matching’’.
  * (76) F. Maltoni, Z. Sullivan and S. Willenbrock, Phys. Rev. D 67 (2003), 093005 doi:10.1103/PhysRevD.67.093005 [arXiv:hep-ph/0301033 [hep-ph]].
  * (77) E. Boos and T. Plehn, Phys. Rev. D 69 (2004), 094005 doi:10.1103/PhysRevD.69.094005 [arXiv:hep-ph/0304034 [hep-ph]].
  * (78) F. Maltoni, T. McElmurry and S. Willenbrock, Phys. Rev. D 72 (2005), 074024 doi:10.1103/PhysRevD.72.074024 [arXiv:hep-ph/0505014 [hep-ph]].
  * (79) F. Maltoni, G. Ridolfi and M. Ubiali, JHEP 07 (2012), 022 [erratum: JHEP 04 (2013), 095] doi:10.1007/JHEP04(2013)095 [arXiv:1203.6393 [hep-ph]].
  * (80) R. V. Harlander, Eur. Phys. J. C 76 (2016) no.5, 252 doi:10.1140/epjc/s10052-016-4093-x [arXiv:1512.04901 [hep-ph]].
  * (81) S. Dittmaier, M. Krämer, A. Muck and T. Schluter, JHEP 03 (2007), 114 doi:10.1088/1126-6708/2007/03/114 [arXiv:hep-ph/0611353 [hep-ph]].
  * (82) S. Dawson, C. B. Jackson and P. Jaiswal, Phys. Rev. D 83 (2011), 115007 doi:10.1103/PhysRevD.83.115007 [arXiv:1104.1631 [hep-ph]].
  * (83) S. Dittmaier, P. Häfliger, M. Krämer, M. Spira and M. Walser, Phys. Rev. D 90 (2014) no.3, 035010 doi:10.1103/PhysRevD.90.035010 [arXiv:1406.5307 [hep-ph]].
  * (84) E. Bagnaschi, H. Bahl, E. Fuchs, T. Hahn, S. Heinemeyer, S. Liebler, S. Patel, P. Slavich, T. Stefaniak and C. E. M. Wagner, et al. Eur. Phys. J. C 79 (2019) no.7, 617 doi:10.1140/epjc/s10052-019-7114-8 [arXiv:1808.07542 [hep-ph]].