# Probing double-aligned two Higgs doublet models at LHC

Shinya Kanemura<EMAIL_ADDRESS>Michihisa Takeuchi
<EMAIL_ADDRESS>Kei Yagyu
<EMAIL_ADDRESS>Department of Physics, Osaka University,
Toyonaka, Osaka 560-0043, Japan

###### Abstract

We consider two Higgs doublet models (THDMs) with both the Higgs potential and
Yukawa interactions being aligned, which we call “double-aligned THDMs”. In
this scenario, coupling constants of the discovered Higgs boson to the
Standard Model (SM) particles are identical to those of the SM Higgs boson,
and flavor changing neutral currents via neutral Higgs bosons do not appear at
tree level. We investigate current constraints and future prospects of the
model by using measurements from flavor experiments and data of multi-lepton
final states at LHC. Especially, we focus on the electroweak pair production
of the additional Higgs bosons with their masses below $2m_{t}$. We find that
the most of the parameter space are already excluded by the current LHC data
when the leptonic decays of the additional Higgs bosons are dominant, which
can be interpreted to the scenario in the Type-X THDM as a special case. We
also clarify the parameter region where the high-luminosity LHC can explore,
and demonstrate the reconstruction of the masses of additional Higgs bosons
from the $b\bar{b}\tau^{+}\tau^{-}$ final states in a few benchmark points.

††preprint: OU-HET-1116

##  I Introduction

The nature of the electroweak (EW) symmetry breaking is one of the biggest
mysteries in the elementary particle physics even after the discovery of the
Higgs boson at LHC. Although the Higgs mechanism in the standard model (SM)
can parameterize the EW symmetry breaking and can explain currently measured
properties of the Higgs boson, the true structure of the Higgs sector is still
unknown as it is not determined by a fundamental principle. In fact, various
possibilities of non-minimal structures of the Higgs sector can be considered
without any contradiction with the current experimental data. Therefore,
comprehensive studies of extended Higgs sectors are important to understand
the nature of EW symmetry breaking as a bottom-up approach.

Two Higgs doublet models (THDMs) are one of the simplest examples of the
extended Higgs sector, which naturally satisfy the electroweak rho parameter
to be unity at tree level. In addition, THDMs are well motivated, because they
appear in various new physics scenarios such as supersymmetry Haber and Kane
(1985), a neutrino-philic model Ma (2001), models with radiative generations
of neutrino masses Zee (1980); Ma (2006); Aoki  _et al._ (2009a), those with
dark matter Barbieri  _et al._ (2006), those with new sources of CP violation
(CPV) and those with strongly first order EW phase transitions Kanemura  _et
al._ (2005); Fromme  _et al._ (2006); Funakubo  _et al._ (1994). In
particular, such CPV and first order EW phase transitions can be used in the
successful scenarios of the electroweak baryogenesis Turok and Zadrozny
(1991); Funakubo  _et al._ (1994); Cline  _et al._ (1996); Fromme  _et al._
(2006); Modak and Senaha (2019); Basler  _et al._ (2021). Thus, among the many
possibilities of non-minimal Higgs sectors, THDMs should phenomenologically be
examined with a special attention.

It has been well known that THDMs generally give rise to flavor changing
neutral currents (FCNCs) via neutral Higgs boson exchanges at tree level. In
order to avoid such FCNCs, a natural flavor conservation scenario Glashow and
Weinberg (1977) has often been considered, where only one of the Higgs
doublets couples to each type of charged fermions, i.e., up-type quarks, down-
type quarks and charged leptons. This scenario can naturally be realized by
introducing a discrete $Z_{2}$ symmetry to the Higgs sector, in which four
independent type of Yukawa interactions appear Grossman (1994); Aoki  _et al._
(2009b); Barger  _et al._ (1990), the so-called Type-I, Type-II, Type-X and
Type-Y. As an alternative possibility to avoid tree level FCNCs, the Yukawa
alignment scenario has also been discussed Pich and Tuzon (2009), in which two
Yukawa matrices for each type of charged fermions are assumed to be
proportional with each other. In this case, Yukawa couplings can generally be
complex, so that we can consider two independent CPV sources in the Yukawa
interactions and the Higgs potential.

Apart from the Yukawa alignment, current LHC data show that measured
properties of the discovered Higgs boson agree with those of the SM Higgs
boson within the uncertainty. This suggests that one of the neutral Higgs
bosons in the THDM, denoting $h$, must have SM-Higgs like properties. There
are two ways to realize such a situation, i.e., taking decoupling or alignment
regions of the parameter space. The former corresponds to the case with all
the masses of additional Higgs bosons being much higher than the EW scale, in
which properties of $h$ approach to those of the SM Higgs boson according to
the decoupling theorem Appelquist and Carazzone (1975); Gunion and Haber
(2003). On the other hand, the latter can be realized by taking the $h$ state
being aligned to be one of the states in the Higgs doublet which possesses the
vacuum expectation value (VEV) same as in the SM 111Such a Higgs doublet can
be defined by taking the so-called Higgs basis Davidson and Haber (2005)
without loss of generality. . Unlike the decoupling case, the scenario with
the “Higgs alignment” provides phenomenologically interesting consequences,
because the additional Higgs bosons can be at the EW scale, which is often
referred as “alignment without decoupling” Aiko  _et al._ (2021). 222In this
case, non-decoupling effects of additional Higgs bosons can be significant,
which provide a sizable deviation in the triple Higgs boson coupling, see
e.g., Refs. Kanemura  _et al._ (2004, 2015); Kanemura and Nagai (2021);
Braathen and Kanemura (2019, 2020). Such a large deviation can be probed in
the double Higgs boson production at future colliders Plehn  _et al._ (1996);
Baur  _et al._ (2002, 2003); Asakawa  _et al._ (2010); Barr  _et al._ (2015);
Papaefstathiou (2015); Cao  _et al._ (2017); Gonçalves  _et al._ (2018);
Biekötter  _et al._ (2019); Abada  _et al._ (2019a, b); Aihara  _et al._
(2019); Klamka (2021); Han  _et al._ (2021); Chiesa  _et al._ (2020).  From
the above mentioned reasons, THDMs with the Yukawa and Higgs alignment provide
a phenomenologically interesting scenario. We call this type of scenario as
“double-alignment”. Recently in Ref. Kanemura  _et al._ (2020), such a double-
aligned scenario has been studied, and found that sizable CPV phases can be
taken without contradiction to bounds from current data of electric dipole
moments (EDMs) Andreev  _et al._ (2018); Abel  _et al._ (2020) as well as
flavor experiments and LHC experiments. See also Refs. Cheung  _et al._
(2020); Abe  _et al._ (2014); Jung and Pich (2014); Altmannshofer  _et al._
(2020); Low  _et al._ (2020) for discussions of EDMs in THDMs. In addition, it
has been shown in Ref. Enomoto  _et al._ (2021) that these CPV phases can be
used to realize a successful scenario of the EW baryogenesis. Furthermore,
effects of CPV can directly be detected by measuring azimuthal angle
distributions in the decay product of additional Higgs bosons Kanemura  _et
al._ (2021) at future lepton colliders such as the International Linear
Collider (ILC) Baer  _et al._ (2013); Asai  _et al._ (2017); Fujii  _et al._
(2017), the Circular Electron Positron Collider (CEPC) Group (2015) and the
Future Circular Collider (FCC-ee) Bicer  _et al._ (2014).

In this paper, we study the possibility of direct searches for additional
Higgs bosons at LHC in the THDM with the double-alignment. In particular, we
focus on the EW pair production of the additional Higgs bosons whose cross
sections are simply determined by the masses of the additional Higgs bosons
Kanemura and Yuan (2002); Cao  _et al._ (2004); Belyaev  _et al._ (2008).
333In Ref. Eberhardt  _et al._ (2021), global fits in the Yukawa aligned THDM
have been performed, where various flavor data, signal strengths of the
discovered Higgs boson and single productions of the additional Higgs bosons
have been combined.  We apply the data for multi-lepton final states at the
LHC Run-II experiment to constrain the parameter space of the THDM, and
combine constraints from flavor experiments such as $B\to X_{s}\gamma$,
$B_{s}\to\mu\mu$ and leptonic tau decays. Furthermore, we extrapolate the Run-
II data to obtain the parameter space expected to be explored at the High-
Luminosity LHC (HL-LHC). We find that a large portion of the parameter space
has already been excluded by the current LHC data, especially the region being
able to be regarded as the Type-X (or lepton specific) THDM.

This paper is organized as follows. In Sec. II, we give the Lagrangian of our
THDM, and define the double-alignment, i.e., the Yukawa alignment and the
Higgs alignment. In Sec. V, we discuss the decays of the additional Higgs
bosons. Sec. IV summarizes various constraints from flavor experiments. Sec. V
is devoted to show our main results, i.e., parameter regions excluded by the
current LHC data and expected to be explored at the HL-LHC. We also
demonstrate the reconstruction of masses of additional Higgs bosons from the
$b\bar{b}\tau^{+}\tau^{-}$ state. Conclusions and discussions are given in
Sec. VI. In Appendix, we show the constraint on the parameter space for the
special case with $\zeta_{u}=0$.

##  II Model

We consider a model with two isospin Higgs doublets without introducing any
symmetry other than the SM gauge symmetry. In this case, the Lagrangian is
invariant under the global $U(2)$ transformation between the two doublets by
using an appropriate redefinition of parameters in the potential and Yukawa
interactions. Using this invariance, we can write down the Lagrangian in terms
of the Higgs basis Davidson and Haber (2005) $\Phi$ and $\Phi^{\prime}$
without loss of generality, which is defined as

$\displaystyle\Phi=\begin{pmatrix}G^{+}\\\
\frac{1}{\sqrt{2}}(v+h^{0}_{1}+iG^{0})\end{pmatrix},\quad\Phi^{\prime}=\begin{pmatrix}H^{+}\\\
\frac{1}{\sqrt{2}}(h^{0}_{2}+ih^{0}_{3})\end{pmatrix},$ (1)

where $v$ is the VEV related to the Fermi constant $G_{F}$ through
$v=(\sqrt{2}G_{F})^{-1/2}$. In Eq. (1), $G^{\pm}$ and $G^{0}$ are the Nambu-
Goldstone bosons which are absorbed into the longitudinal component of the
$W^{\pm}$ and $Z$ boson, respectively, while $H^{\pm}$ and
$h_{i}^{0}~{}(i=1,2,3)$ are the physical charged and neutral Higgs bosons,
respectively.

The potential is written in the Higgs basis as

$\displaystyle V=$ $\displaystyle
m^{2}|\Phi|^{2}+M^{2}|\Phi^{\prime}|^{2}-\left(\mu^{2}\Phi^{\dagger}\Phi^{\prime}+\text{h.c.}\right)+\frac{\lambda_{1}}{2}|\Phi|^{4}+\frac{\lambda_{2}}{2}|\Phi^{\prime}|^{4}+\lambda_{3}|\Phi|^{2}|\Phi^{\prime}|^{2}+\lambda_{4}|\Phi^{\dagger}\Phi^{\prime}|^{2}$
$\displaystyle+\left[\frac{\lambda_{5}}{2}(\Phi^{\dagger}\Phi^{\prime})+\lambda_{6}|\Phi|^{2}+\lambda_{7}|\Phi^{\prime}|^{2}\right](\Phi^{\dagger}\Phi^{\prime})+\text{h.c.},$
(2)

where $\mu^{2}$ and $\lambda_{5,6,7}$ are complex in general. The tadpole
conditions, vanishment of the linear term of $h_{i}^{0}$, provide

$\displaystyle
m^{2}=\frac{1}{2}\lambda_{1}v^{2},\quad\mu^{2}=\frac{1}{2}\lambda_{6}v^{2}.$
(3)

We note that the second equation is given by the conditions with respect to
$h_{2}^{0}$ or $h_{3}^{0}$. By imposing the tadpole condition, the squared
mass of $H^{\pm}$ is given by

$\displaystyle m_{H^{\pm}}^{2}=M^{2}+\frac{1}{2}\lambda_{3}v^{2}.$ (4)

The squared-mass matrix for the neutral Higgs bosons in the basis of
$(h^{0}_{1},h^{0}_{2},h^{0}_{3})$ is given by

$\displaystyle\mathcal{M}^{2}=v^{2}\begin{pmatrix}\lambda_{1}&\textrm{Re}[\lambda_{6}]&-\textrm{Im}[\lambda_{6}]\\\
\textrm{Re}[\lambda_{6}]&\frac{M^{2}}{v^{2}}+\frac{1}{2}(\lambda_{3}+\lambda_{4}+\textrm{Re}[\lambda_{5}])&-\frac{1}{2}\textrm{Im}[\lambda_{5}]\\\
-\textrm{Im}[\lambda_{6}]&-\frac{1}{2}\textrm{Im}[\lambda_{5}]&\frac{M^{2}}{v^{2}}+\frac{1}{2}(\lambda_{3}+\lambda_{4}-\textrm{Re}[\lambda_{5}])\end{pmatrix}.$
(5)

The mass eigenstates $H_{i}^{0}~{}(i=1,2,3)$ can be defined by introducing the
orthogonal matrix $R$ as $h_{i}^{0}=O_{ij}H_{j}^{0}$ with the eigenvalues
$\text{diag}(m_{H_{1}^{0}}^{2},m_{H_{2}^{0}}^{2},m_{H_{3}^{0}}^{2})\equiv
O^{T}{\cal M}^{2}O$. We define $m_{H_{1}^{0}}\leq m_{H_{2}^{0}}\leq
m_{H_{3}^{0}}^{2}$, and identify the $H_{1}^{0}$ state with the discovered
Higgs boson with the mass of 125 GeV.

In this paper, we impose the alignment condition:

$\displaystyle\lambda_{6}=0,$ (6)

by which the $h_{1}^{0}$ state coincides with the mass eigenstate $H_{1}^{0}$,
and as a result couplings of $H_{1}^{0}$ to the gauge bosons and fermions
agree with those of the SM Higgs boson at tree level. We refer this condition
as “Higgs alignment”. By rephasing $\Phi^{\prime}$, the complex phase of
$\lambda_{5}$ can be removed without loss of generality, so that the mass
matrix becomes diagonal form. We note that at this stage we cannot identify
the additional Higgs boson $H_{2}^{0}$ ($H_{3}^{0}$) with the CP-even (CP-odd)
state, because its CP property depends on the structure of the Yukawa
interaction.

In this scenario, there are 7 free parameters444The number of parameters can
consistently be counted as follows: 14 (initial number of real parameters in
the potential) $-$ 3 (tadpole conditions) $-$ 2 (Higgs alignment) $-$ 1
(rephasing) = 7 [Eq. (7)] + 1 ($m_{H_{1}^{0}}=125$ GeV).  which can be chosen
as follows

$\displaystyle
M^{2},~{}m_{H^{0}_{2}},~{}m_{H^{0}_{3}},~{}m_{H^{\pm}},~{}\lambda_{2},~{}|\lambda_{7}|~{}\textrm{and}~{}\theta_{7},$
(7)

where $\theta_{7}\equiv\arg[\lambda_{7}]\in(-\pi,\pi]$.

The most general Yukawa interactions are given in the mass eigenstates of
fermions as

$\displaystyle\mathcal{L}_{\rm Y}=$
$\displaystyle-\bar{Q}_{L}^{u}\left(\sqrt{2}\frac{M_{u}}{v}\tilde{\Phi}+\rho_{u}\tilde{\Phi}^{\prime}\right)u_{R}-\bar{Q}_{L}^{d}\left(\sqrt{2}\frac{M_{d}}{v}\Phi+\rho_{d}\Phi^{\prime}\right)d_{R}$
$\displaystyle-\bar{L}_{L}\left(\sqrt{2}\frac{M_{e}}{v}\Phi+\rho_{e}\Phi^{\prime}\right)e_{R}+\text{h.c.},$
(8)

where $Q_{L}^{u}=(u_{L},V_{\rm CKM}d_{L})^{T}$, $Q_{L}^{d}=(V_{\rm
CKM}^{\dagger}u_{L},d_{L})^{T}$ and $L_{L}=(\nu_{L},e_{L})$ are left-handed
doublet fermions with $V_{\rm CKM}$ being the Cabibbo-Kobayashi-Maskawa (CKM)
matrix, while $u_{R}$, $d_{R}$ and $e_{R}$ are right-handed up-type quarks,
down-type quarks and charged leptons, respectively. The charge conjugation of
the Higgs doublets is denoted as
$\tilde{\Phi}^{(\prime)}=i\tau_{2}\Phi^{(\prime)*}$. In the first term of each
parentheses, $M_{f}$ ($f=u,d,e$) denote diagonalized mass matrices, while in
the second term $\rho_{f}$ represent arbitrary $3\times 3$ complex matrices.
In this expression, we do not explicitly show flavor indices.

It is clear that the $\rho_{f}$ terms give rise to tree-level FCNCs mediated
by neutral Higgs bosons, particularly via $H_{2}^{0}$ and/or $H_{3}^{0}$ in
the Higgs alignment limit defined in Eq. (6). In order to avoid such FCNCs, we
impose so-called the Yukawa alignment Pich and Tuzon (2009), i.e.,

$\displaystyle\rho_{f}=\sqrt{2}\zeta_{f}\frac{M_{f}}{v},$ (9)

where $\zeta_{f}$ are arbitrary complex parameters. The Yukawa alignment can
also be described in the general basis of the two doublets as the assumption
that two Yukawa matrices for a fermion type $f$ are proportional to each
other. We note that the $\zeta_{f}$ parameters are flavor universal, e.g.,
$\zeta_{\mu}=\zeta_{\tau}=\zeta_{e}$, due to the above assumption. Because the
$\zeta_{f}$ parameters are complex, new sources of the CPV appear in addition
to the potential parameter $\lambda_{7}$. In Ref. Kanemura  _et al._ (2020),
it has been shown that we can take sizable CPV phases, while severe
constraints from experiments of the electron EDM Andreev  _et al._ (2018) and
the neutron EDM Abel  _et al._ (2020) can be avoided by using cancellation
among Barr-Zee diagrams Barr and Zee (1990) with fermion and scalar boson
loops due to the independent phases from Yukawa interactions and the Higgs
potential. Such a scenario is compatible with constraints from the other
flavor experiments as well as the current LHC data. We note that in THDMs with
a softly-broken $Z_{2}$ symmetry Glashow and Weinberg (1977) these $\zeta_{f}$
parameters are determined by one parameter $\tan\beta$ (the ratio of the two
Higgs VEVs) depending on the type of Yukawa interactions Barger  _et al._
(1990); Aoki  _et al._ (2009b) by

$\displaystyle\begin{split}(\zeta_{u},\zeta_{d},\zeta_{e})&=(\cot\beta,\cot\beta,\cot\beta)\quad\quad\quad\text{in
Type-I},\\\
(\zeta_{u},\zeta_{d},\zeta_{e})&=(\cot\beta,-\tan\beta,-\tan\beta)\quad\text{in
Type-II},\\\
(\zeta_{u},\zeta_{d},\zeta_{e})&=(\cot\beta,\cot\beta,-\tan\beta)\quad\quad\text{in
Type-X},\\\
(\zeta_{u},\zeta_{d},\zeta_{e})&=(\cot\beta,-\tan\beta,\cot\beta)\quad\quad\text{in
Type-Y}.\end{split}$ (10)

In the Higgs and Yukawa alignment, the double-alignment, the Higgs boson
couplings are expressed as follows

$\displaystyle{\cal L}_{\rm int}$
$\displaystyle=-\sum_{f=u,d,e}\bar{f}\frac{M_{f}}{v}\left[H_{1}^{0}+|\zeta_{f}|(\cos\theta_{f}+i\sin\theta_{f}\gamma_{5})H_{2}^{0}+2I_{f}|\zeta_{f}|(\sin\theta_{f}-i\cos\theta_{f}\gamma_{5})H_{3}^{0}\right]f$
$\displaystyle-\frac{\sqrt{2}}{v}\bar{u}\left(P_{R}\zeta_{d}V_{\rm
CKM}M_{d}-P_{L}\zeta_{u}^{*}M_{u}^{\dagger}V_{\rm
CKM}^{\dagger}\right)dH^{+}-\frac{\sqrt{2}}{v}\bar{\nu}\left(P_{R}\zeta_{e}M_{e}\right)eH^{+}+\text{h.c.},$
(11)

where $\theta_{f}\equiv\arg[\zeta_{f}]\in(-\pi,\pi]$, $I_{u}=1/2$,
$I_{d}=I_{e}=-1/2$ and $P_{L}$ ($P_{R}$) is the projection operator for left-
(right) handed fermions. It is clear that non-zero phases $\theta_{f}$ result
in the CP mixing in the Yukawa sector. We note that the $H_{1}^{0}VV$
($V=W,Z$) couplings coincide with the SM values at tree level due to the Higgs
alignment while $H_{2,3}^{0}VV$ couplings vanish. On the other hand, there are
Higgs-Higgs-Gauge type interactions for the additional Higgs bosons such as
$H_{2}^{0}H_{3}^{0}Z$ and $H_{2,3}^{0}H^{\pm}W^{\mp}$, which is
phenomenologically important for the EW pair productions of the additional
Higgs bosons and their decays as discussed in the following sections.

##  III Decays of the additional Higgs bosons

We discuss decays of additional Higgs bosons in the double-alignment limit,
where they can mainly decay into a fermion pair or a lighter additional Higgs
boson with a gauge boson as long as these are kinematically allowed. In the
following discussion, we focus on the case where the charged Higgs boson mass
is degenerate with one of the additional neutral Higgs bosons in order to
avoid the $T$ parameter constraint Peskin and Takeuchi (1990); Bertolini
(1986); Peskin and Wells (2001); Grimus  _et al._ (2008); Kanemura  _et al._
(2011). In this case, there are two possible scenarios denoted as the light
$H^{\pm}$ scenario: $m_{H^{\pm}}=m_{H_{2}^{0}}$ and the heavy $H^{\pm}$
scenario: $m_{H^{\pm}}=m_{H_{3}^{0}}$, where we define $m_{H_{3}^{0}}\geq
m_{H_{2}^{0}}\geq m_{H_{1}}$ as mentioned in Sec. II.

The decay rates of the Higgs bosons into a fermion pair are given by

$\displaystyle\Gamma(H_{2,3}^{0}\to f\bar{f})$
$\displaystyle=\frac{\sqrt{2}G_{F}m_{H_{2,3}}^{3}}{8\pi}N_{f}^{c}x_{f}|\zeta_{f}|^{2}[1-2x_{f}(1\pm\cos
2\theta_{f})]\lambda^{1/2}\left(x_{f},x_{f}\right),$ (12)
$\displaystyle\Gamma(H^{\pm}\to ff^{\prime})$
$\displaystyle=\frac{\sqrt{2}G_{F}m_{H^{\pm}}^{3}}{8\pi}|V_{\rm
CKM}^{ff^{\prime}}|^{2}N_{f}^{c}\lambda^{1/2}(x_{f},x_{f^{\prime}})$
$\displaystyle\times\left[(x_{f}|\zeta_{f}|^{2}+x_{f^{\prime}}|\zeta_{f^{\prime}}|^{2})(1-x_{f}-x_{f^{\prime}})+4x_{f}x_{f^{\prime}}\text{Re}(\zeta_{f}^{*}\zeta_{f^{\prime}})\right],$
(13)

where $\lambda(x,y)=(1-x-y)^{2}-4xy$, $x_{i}=m_{i}^{2}/m_{\cal H}^{2}$ with
$m_{\cal H}$ being the mass of the decaying Higgs boson and $N_{f}^{c}=3(1)$
for $f$ being quarks (leptons). As we mentioned in Sec. II, $\zeta_{f}$
factors are assumed to be flavor universal, e.g.,
$\zeta_{\mu}=\zeta_{\tau}=\zeta_{e}$. For the case with non-zero mass
differences among the additional Higgs bosons, the following decay rates have
to be added to their total widths:

$\displaystyle\Gamma(H^{0}_{3}\to\phi
V^{(\ast)})=\begin{cases}\frac{\sqrt{2}G_{F}m_{H_{3}^{0}}^{3}}{16\pi}\,\lambda^{3/2}(x_{\phi},x_{V})\quad(m_{H_{3}^{0}}-m_{\phi}\geq
m_{V})\\\
\frac{9m_{H_{3}^{0}}}{16\pi^{3}}G_{F}^{2}m_{V}^{4}\delta_{V}\,G(x_{\phi},x_{V})\quad(m_{H_{3}^{0}}-m_{\phi}<m_{V})\end{cases},$
(14) $\displaystyle\Gamma(H^{\pm}\to
H_{2}^{0}W^{\pm(\ast)})=\begin{cases}\frac{\sqrt{2}G_{F}m_{H^{\pm}}^{3}}{16\pi}\,\lambda^{3/2}(x_{H^{\pm}},x_{W})\quad(m_{H^{\pm}}-m_{H_{2}^{0}}\geq
m_{W})\\\
\frac{9m_{H^{\pm}}}{16\pi^{3}}G_{F}^{2}m_{W}^{4}G(x_{H_{2}},x_{W})\quad(m_{H^{\pm}}-m_{H_{2}^{0}}<m_{W})\end{cases},$
(15)

where 555In Eq. (14) with $\phi=H^{\pm}$ and $V=W$, this expression does not
take the sum over two possible states, i.e., $H^{+}W^{-(*)}$ and
$H^{-}W^{+(*)}$.  $\phi=H^{0}_{2}\,(H^{\pm})$ for
$V^{(*)}=Z^{(*)}\,(W^{\mp(*)})$, $\delta_{W}=1$ and
$\delta_{Z}=\frac{7}{6}-\frac{20}{9}\sin^{2}\theta_{W}+\frac{80}{27}\sin^{4}\theta_{W}$.
The function $G$ is the phase space function for three body decays, see e.g.,
Ref. Aoki  _et al._ (2012) for the explicit form of this function:

$\displaystyle G(x,y)$
$\displaystyle=\frac{1}{12y}\Bigg{\\{}2\left(x-1\right)^{3}-9\left(x^{2}-1\right)y+6\left(x-1\right)y^{2}-3\left[1+\left(x-y\right)^{2}-2y\right]y\log
x$
$\displaystyle+6\left(1+x-y\right)y\sqrt{-\lambda(x,y)}\left[\arctan\left(\frac{-1+x-y}{\sqrt{-\lambda(x,y)}}\right)+\arctan\left(\frac{-1+x+y}{\sqrt{-\lambda(x,y)}}\right)\right]\Bigg{\\}}.$
(16)

We note that loop induced decays of the additional neutral Higgs bosons, i.e.,
$H_{2,3}^{0}\to gg/\gamma\gamma/Z\gamma$ are also taken into account in
numerical evaluations, while those of the charged Higgs boson, i.e.,
$H^{\pm}\to W^{\pm}\gamma$ Capdequi Peyranere  _et al._ (1991) and $H^{\pm}\to
W^{\pm}Z$ Kanemura (2000) are neglected because of their tiny partial widths.

As we mentioned in the above, the decay modes of $H_{3}^{0}$ can be classified
into two categories, i.e., fermionic modes ($H_{3}^{0}\to f\bar{f}$) and
bosonic modes ($H_{3}^{0}\to\phi V^{(*)}$). The relative size of the branching
ratios for the fermionic and bosonic modes can be expressed by introducing the
following ratio $R$:

$\displaystyle R$ $\displaystyle\equiv\frac{\sum_{f}\Gamma(H_{3}^{0}\to
f\bar{f})}{\sum_{f}\Gamma(H_{3}^{0}\to
f\bar{f})+\sum_{V}\Gamma(H_{3}^{0}\to\phi V^{(*)})}\simeq\sum_{f}{\cal
B}(H_{3}^{0}\to f\bar{f}).$ (17)

In addition, we introduce

$\displaystyle R_{\tau}$
$\displaystyle\equiv\frac{\Gamma(H_{3}^{0}\to\tau^{+}\tau^{-})}{\sum_{f}\Gamma(H_{3}^{0}\to
f\bar{f})},$ (18)

by which we can parameterize the relative magnitude of the branching ratio of
$H_{3}^{0}\to\tau^{+}\tau^{-}$ among the fermionic modes. These $R$ parameters
can simply be rewritten as

$\displaystyle R=\frac{1}{1+r/\zeta^{2}},\quad
R_{\tau}=\frac{|\zeta_{e}|^{2}}{\zeta^{2}},$ (19)

where

$\displaystyle\zeta^{2}$ $\displaystyle=\frac{\sum_{f}\Gamma(H_{3}^{0}\to
f\bar{f})}{\Gamma_{0}},\quad\Gamma_{0}=\frac{\sqrt{2}G_{F}}{8\pi}m_{H_{3}^{0}}m_{\tau}^{2},$
(20) $\displaystyle r$
$\displaystyle=\begin{cases}\frac{m_{H_{3}^{0}}^{2}}{2m_{\tau}^{2}}\sum_{V,\phi}\lambda^{3/2}\left(\frac{m_{\phi}^{2}}{m_{H_{3}^{0}}^{2}},\frac{m_{V}^{2}}{m_{H_{3}^{0}}^{2}}\right)\quad(m_{H_{3}^{0}}-m_{\phi}\geq
m_{V})\\\
\frac{9}{2\sqrt{2}\pi^{2}}\frac{G_{F}}{m_{\tau}^{2}}\sum_{V,\phi}m_{V}^{4}\delta_{V}G\left(\frac{m_{\phi}^{2}}{m_{H_{3}^{0}}^{2}},\frac{m_{V}^{2}}{m_{H_{3}^{0}}^{2}}\right)\quad(m_{H_{3}^{0}}-m_{\phi}<m_{V}),\end{cases}$
(21)

with the summation $\sum_{V,\phi}$ being taken to be $(V,\phi)=(Z,H_{3}^{0})$
and $(W,H^{\pm})$. In particular, for the case with $m_{b}\ll
m_{H_{2,3}}<2m_{t}$, which will be mainly considered in Sec. V, $\zeta^{2}$
takes a significantly simple form as

$\displaystyle\zeta^{2}$ $\displaystyle\simeq\frac{1}{m_{\tau}^{2}}\sum_{f\neq
t}m_{f}^{2}N_{f}^{c}|\zeta_{f}|^{2}.$ (22)

We note that $H_{2}^{0}$ can only decay into a fermion pair at tree level, so
that the $R$ value for $H_{2}^{0}$, i.e., $R|_{H_{3}^{0}\to H_{2}^{0}}$ is
unity.

Similar to the neutral Higgs bosons, we define the ratio parameters for
$H^{\pm}$ as follows:

$\displaystyle R^{\pm}$ $\displaystyle\equiv\frac{\sum_{f}\Gamma(H^{\pm}\to
f\bar{f}^{\prime})}{\sum_{f}\Gamma(H^{\pm}\to
f\bar{f}^{\prime})+\Gamma(H^{\pm}\to H_{2}^{0}W^{\pm(*)})}\simeq\sum_{f}{\cal
B}(H^{\pm}\to f\bar{f}^{\prime}),$ (23) $\displaystyle R_{\tau}^{\pm}$
$\displaystyle\equiv\frac{\Gamma(H^{\pm}\to\tau\nu)}{\sum_{f}\Gamma(H^{\pm}\to
f\bar{f}^{\prime})}.$ (24)

By introducing the quantities $\zeta_{\pm}^{2}$ and $r_{\pm}$, these $R$
parameters can be expressed as

$\displaystyle R^{\pm}=\frac{1}{1+r_{\pm}/\zeta_{\pm}^{2}},\quad
R_{\tau}^{\pm}=\frac{|\zeta_{e}|^{2}}{\zeta_{\pm}^{2}},$ (25)

where

$\displaystyle\zeta_{\pm}^{2}$ $\displaystyle=\frac{\sum_{f}\Gamma(H^{\pm}\to
f\bar{f}^{\prime})}{\Gamma_{0}|_{H_{3}^{0}\to
H^{\pm}}}\simeq|\zeta_{e}|^{2}+3\left(1-\frac{m_{t}^{2}}{m_{H^{\pm}}^{2}}\right)^{2}\left(\frac{m_{t}^{2}}{m_{\tau}^{2}}|\zeta_{u}|^{2}+\frac{m_{b}^{2}}{m_{\tau}^{2}}|\zeta_{d}|^{2}\right),$
(26) $\displaystyle r_{\pm}$
$\displaystyle=\begin{cases}\frac{m_{H^{\pm}}^{2}}{2m_{\tau}^{2}}\lambda^{3/2}\left(\frac{m_{H_{2}^{0}}^{2}}{m_{H^{\pm}}^{2}},\frac{m_{W}^{2}}{m_{H^{\pm}}^{2}}\right)\quad(m_{H^{\pm}}-m_{H_{2}^{0}}\geq
m_{W})\\\
\frac{9}{2\sqrt{2}\pi^{2}}\frac{G_{F}}{m_{\tau}^{2}}m_{W}^{4}G\left(\frac{m_{H_{2}^{0}}^{2}}{m_{H^{\pm}}^{2}},\frac{m_{W}^{2}}{m_{H^{\pm}}^{2}}\right)\quad(m_{H^{\pm}}-m_{H_{2}^{0}}<m_{W})\end{cases},$
(27)

Unlike the $H_{3}^{0}$ decays, only the decay of $H^{\pm}\to
H_{2}^{0}W^{\pm(*)}$ is allowed for the bosonic decay mode.

In terms of these $R$ parameters, the decay branching ratios into a tau lepton
pair, which will be important in the discussion for the phenomenology at LHC,
can simply be expressed as

$\displaystyle{\cal B}(H_{2}^{0}\to\tau^{+}\tau^{-})\simeq R_{\tau},\quad{\cal
B}(H_{3}^{0}\to\tau^{+}\tau^{-})\simeq RR_{\tau},\quad{\cal
B}(H^{\pm}\to\tau^{\pm}\nu)\simeq R^{\pm}R_{\tau}^{\pm}.$ (28)

It would be important to discuss the critical values of $\zeta$ and
$\zeta_{\pm}$ which give the sum of the fermionic decay branching ratios to be
50%. Such critical values, denoted as $\zeta_{50}$ and $\zeta_{\pm 50}$, can
be expressed as $\zeta_{50}=\sqrt{r}$ and $\zeta_{\pm 50}=\sqrt{r_{\pm}}$ from
Eqs. (19) and (25), respectively. Since $r$ and $r_{\pm}$ depend only on the
masses of additional Higgs bosons, $\zeta_{50}$ and $\zeta_{\pm 50}$ are
determined as a function of these masses.

Figure 1: Values of $\zeta_{50}$ (left) and $\zeta_{\pm 50}$ (right) as a
function of the mass difference $m_{H_{3}^{0}}-m_{H_{2}^{0}}$ for
$m_{H_{2}^{0}}=100$ GeV (black), 300 GeV (blue) and 500 GeV (red). For the
left panel, the solid and dashed curves show the case with the heavy
$(m_{H^{\pm}}=m_{H_{3}^{0}})$ and the light ($m_{H^{\pm}}=m_{H_{2}^{0}}$)
$H^{\pm}$ scenario, respectively.

Fig. 1 shows the values of $\zeta_{50}$ (left panel) and $\zeta_{\pm 50}$
(right panel) as a function of $m_{H_{3}^{0}}-m_{H_{2}^{0}}$. We see that for
$m_{H_{3}^{0}}-m_{H_{2}^{0}}\ll m_{W}$ both $\zeta_{50}$ and $\zeta_{\pm 50}$
are much smaller than unity, because the bosonic decay modes are significantly
suppressed by the phase space, while they can be of order one for
$m_{H_{3}^{0}}-m_{H_{2}^{0}}\sim m_{W}$. When the on-shell decays open,
$\zeta_{50}$ and $\zeta_{\pm 50}$ are of order 10 or larger. We also see that
$\zeta_{50}$ and $\zeta_{\pm 50}$ almost do not depend on $m_{H_{3}^{0}}$, so
that they are essentially determined by the mass difference
$m_{H_{3}^{0}}-m_{H_{2}^{0}}$. It is clear that larger values of $\zeta_{50}$
are required for $m_{H^{\pm}}=m_{H_{2}^{0}}$ (dashed curves) as compared with
that with $m_{H^{\pm}}=m_{H_{3}^{0}}$ (solid curves), because the
$H_{3}^{0}\to H^{\pm}W^{\pm(*)}$ mode also contribute to the $r$ value.

Figure 2: Contour plots for $R_{\tau}=50\%$ (left) and $R_{\tau}^{\pm}=50\%$
(left) on the $|\zeta_{d}|/|\zeta_{e}|$ and $|\zeta_{u}|/|\zeta_{e}|$ plane
for each fixed value of $m_{H_{2}^{0}}$ (left) and $m_{H^{\pm}}$ (right).

Similarly, we discuss the critical values of $\zeta_{f}$ which lead to
$R_{\tau}=50\%$ or $R_{\tau}^{\pm}=50\%$. From Eqs. (19) and (25), $R_{\tau}$
and $R_{\tau}^{\pm}$ are determined by the ratios of the $\zeta_{f}$
parameters, i.e., $\zeta_{u}/\zeta_{e}$ and $\zeta_{d}/\zeta_{e}$. More
concretely, the inverse of $R_{\tau}$ and $R_{\tau}^{\pm}$ can be expressed
for $\theta_{f}=0$ as

$\displaystyle R_{\tau}^{-1}$ $\displaystyle\simeq
1+\frac{3m_{b}^{2}}{m_{\tau}^{2}}\frac{|\zeta_{d}|^{2}}{|\zeta_{e}|^{2}}+\left[\frac{3m_{c}^{2}}{m_{\tau}^{2}}+\theta_{tt}\frac{3m_{t}^{2}}{m_{\tau}^{2}}\left(1-\frac{4m_{t}^{2}}{m_{H_{3}^{0}}^{2}}\right)^{3/2}\right]\frac{|\zeta_{u}|^{2}}{|\zeta_{e}|^{2}},$
(29) $\displaystyle(R_{\tau}^{\pm})^{-1}$ $\displaystyle\simeq
1+\left[\frac{3m_{s}^{2}}{m_{\tau}^{2}}+\theta_{tb}\frac{3m_{b}^{2}}{m_{\tau}^{2}}\left(1-\frac{m_{t}^{2}}{m_{H^{\pm}}^{2}}\right)^{2}\right]\frac{|\zeta_{d}|^{2}}{|\zeta_{e}|^{2}}+\left[\frac{3m_{c}^{2}}{m_{\tau}^{2}}+\theta_{tb}\frac{3m_{t}^{2}}{m_{\tau}^{2}}\left(1-\frac{m_{t}^{2}}{m_{H^{\pm}}^{2}}\right)^{2}\right]\frac{|\zeta_{u}|^{2}}{|\zeta_{e}|^{2}},$
(30)

where $\theta_{tt}\equiv\theta(m_{H_{3}^{0}}-2m_{t})$ and
$\theta_{tb}\equiv\theta(m_{H^{\pm}}-m_{t}-m_{b})$ with $\theta(x)$ being the
Heaviside step function. These expressions show that the values of
$(R_{\tau})^{-1}-1$ and $(R_{\tau}^{\pm})^{-1}-1$ obey the equation of an
ellipse with the $x$ ($y$) axis corresponding to $|\zeta_{d}|/|\zeta_{e}|$
$(|\zeta_{u}|/|\zeta_{e}|)$. This is numerically shown in Fig. 2. Here, we
show the first quadrant of the ellipse. The points on each curve represent the
required values of $|\zeta_{d}|/|\zeta_{e}|$ and $|\zeta_{u}|/|\zeta_{e}|$ to
obtain $R_{\tau}=50\%$ (left panel) and $R_{\tau}^{\pm}=50\%$ (right panel).
For the case below the top threshold, the curve almost does not depend on
$m_{H_{2}^{0}}$, and the interceptions of $x$ and $y$ axis are simply
determined by the fermion mass ratio extracted from Eqs. (29) and (30), i.e.,
$|\zeta_{d}|/|\zeta_{e}|\simeq m_{\tau}/(\sqrt{3}m_{b})\simeq 0.35$ and
$|\zeta_{u}|/|\zeta_{e}|\simeq m_{\tau}/(\sqrt{3}m_{c})\simeq 1.6$. On the
other hand, above the top threshold, the dependence of $m_{H_{2}^{0}}$
slightly appears due to the phase function as seen in Eqs. (29) and (30). For
$R_{\tau}^{\pm}$ (right panel), both the $|\zeta_{d}|/|\zeta_{e}|$ and
$|\zeta_{u}|/|\zeta_{e}|$ depend on $m_{H_{2}^{0}}$, because the phase space
suppression is multiplied to both the $|\zeta_{d}|^{2}/|\zeta_{e}|^{2}$ and
$|\zeta_{u}|^{2}/|\zeta_{e}|^{2}$ terms for $m_{H^{\pm}}>m_{t}+m_{b}$.

By looking at Figs. 1 and 2 and using Eq. (28), the Higgs boson decays into a
tau lepton pair can well be estimated.

##  IV Flavor Constraints

$m_{t}$ [GeV] |  $m_{b}$ [GeV] |  $m_{c}$ [GeV] |  $m_{Z}$ [GeV] |  $m_{W}$ [GeV] |  $m_{H_{1}}$ [GeV] | $\alpha_{s}(m_{Z})$ | $\alpha_{\rm em}^{-1}$ |
---|---|---|---|---|---|---|---|---
172.76 Zyla  _et al._ (2020) | 4.561 Alberti  _et al._ (2015) | 1.092 Alberti  _et al._ (2015) | 91.1876 Zyla  _et al._ (2020) | 80.379 Zyla  _et al._ (2020) | 125.25 Zyla  _et al._ (2020) | 0.1182 Aoki  _et al._ (2020) | 137.036 Zyla  _et al._ (2020) |
${\cal B}(\bar{B}\to X_{c}\ell^{-}\nu)$ | $|V_{ts}^{*}V_{tb}/V_{cb}|^{2}$ | $|V_{cb}|$ |  $\bar{\tau}(B_{s}^{0})$ [ps] |  $\tau(B_{sL}^{0})$ [ps] |  $\tau(B_{sH}^{0})$ [ps] |  $\Delta\Gamma_{s}$ [ps-1] |  $m_{B_{s}}$ [GeV] |  $f_{B_{s}}$ [GeV]
0.1065 Amhis  _et al._ (2021) | 0.9626 Charles  _et al._ (2015); Czakon  _et al._ (2015) | 0.041 Charles  _et al._ (2015) | 1.510 Amhis  _et al._ (2021) | 1.414 Amhis  _et al._ (2021) | 1.619 Amhis  _et al._ (2021) | 0.090 Amhis  _et al._ (2021) | 5.367 Zyla  _et al._ (2020) | 0.2284 Aoki  _et al._ (2020)
Table 1: Input SM parameters. We take the central values of these parameters.

In our scenario, the additional Higgs bosons can mediate various flavor
processes and change their predictions from those in the SM. Such effects can
be translated into constraints on the masses of Higgs bosons and the
$\zeta_{f}$ parameters, defined in Sec. II.

Let us first discuss the constraint from $B\to X_{s}\gamma$ decay process. The
current experimental value of the branching ratio is given by Amhis  _et al._
(2021)

$\displaystyle{\cal B}(B\to X_{s}\gamma)=(3.32\pm 0.15)\times 10^{-4}.$ (31)

In THDMs, the charged Higgs boson can run in the one-loop diagram instead of
the W boson, so that there is a sensitivity to $m_{H^{\pm}}$ and quark Yukawa
couplings which correspond to $\zeta_{u}$ and $\zeta_{d}$ in our scenario. It
has been known that QCD corrections can sizably change the prediction, e.g.,
about a few 10% at NLO with respect to the LO prediction in the SM Ciuchini
_et al._ (1998). The NNLO QCD corrections have also been calculated in Ref.
Misiak and Steinhauser (2017) in the THDM. We implement the NLO QCD and QED
corrections to the decay rate according to Ref. Borzumati and Greub (1998) and
Kagan and Neubert (1999), respectively. 666We set the renormalization scale of
the $B$ meson decay $\mu_{b}$ to be the pole mass of the bottom quark $m_{b}$,
while that of the matching to the full EW theory to be $m_{W}$.  The
constraint from $B\to X_{s}\gamma$ is particularly important for larger values
of $\zeta_{u}$ which correspond to the case with a smaller $\tan\beta$ value
in the THDMs with the softly-broken $Z_{2}$ symmetry, see e.g., Misiak and
Steinhauser (2017).

In addition to $B\to X_{s}\gamma$, we take into account the constraint from
$B_{s}\to\mu\mu$ decay. The current experimental value of the branching ratio
is given by Amhis  _et al._ (2021)

$\displaystyle{\cal B}(B_{s}\to\mu\mu)=(3.1\pm 0.6)\times 10^{-9}.$ (32)

In the aligned THDM, the charged Higgs boson can contribute to the process via
box and penguin type diagrams. The latter also contains neutral Higgs
exchanges. Unlike the $B\to X_{s}\gamma$ process, the $B_{s}\to\mu\mu$ process
has a sensitivity to the masses of neutral Higgs bosons and the $\zeta_{e}$
parameters in addition to $m_{H^{\pm}}$, $\zeta_{u}$ and $\zeta_{d}$.
Therefore, even in a case with smaller $\zeta_{u}$ which is not excluded by
$B\to X_{s}\gamma$, scenarios with larger $\zeta_{d}$ and/or $\zeta_{e}$ can
be excluded. We implement the decay rate of $B_{s}\to\mu\mu$ at one-loop level
according to Ref. Li  _et al._ (2014) by utilizing the Wilson coefficient
$C_{10}(\bar{s}\gamma_{\mu}P_{L}b)(\bar{\mu}\gamma^{\mu}\gamma_{5}\mu)$ in the
SM evaluated at NNLO of QCD as Bobeth  _et al._ (2014)

$\displaystyle C_{10}^{\rm
SM}=-0.938\times\left(\frac{m_{t}}{173.1~{}\text{GeV}}\right)^{1.53}\times\left(\frac{\alpha_{s}(m_{Z})}{0.1184}\right)^{-0.09}.$
(33)

Finally, we take into account the leptonic tau decays which are mediated by
the charged Higgs boson in addition to the contribution from the $W$ boson at
tree level. Thus, the combination of $\zeta_{e}^{2}/m_{H^{\pm}}^{2}$ is
constrained Krawczyk and Temes (2005). Because the coupling of $H^{\pm}$ with
an electron is negligibly small, the contributions of $H^{\pm}$ to $\tau\to
e\nu_{\tau}\bar{\nu}_{e}$ and $\mu\to e\nu_{\mu}\bar{\nu}_{e}$ are negligible.
This causes the violation of lepton flavor universality, which has been
stringently constrained by experiments. The current world average is Amhis
_et al._ (2021)

$\displaystyle(g_{\tau}/g_{\mu})_{\rm ave}=0.9999\pm 0.0014.$ (34)

The left-hand side of the above expression can be compared with the theory
prediction of $\Gamma(\tau\to\mu\nu_{\tau}\bar{\nu}_{\mu})_{\rm
THDM}/\Gamma(\tau\to\mu\nu_{\tau}\bar{\nu}_{\mu})_{\rm SM}$.

For numerical evaluations, we use the SM input parameters shown in Table 1,
and show the parameter region constrained by the flavor data in the next
section.

##  V Direct searches for additional Higgs bosons at LHC

We discuss constraints on the model parameters from direct searches for the
additional Higgs bosons by using the current LHC data in the double-aligned
THDM. We also study the parameter space expected to be explored at the HL-LHC.
For simplicity, we here consider the CP-conserving case, i.e., $\theta_{f}=0$,
in which the constraint from LHC data almost does not depend on these phases
while that from flavor experiments discussed in the previous section has a
sensitivity to these phases. In addition, we concentrate on the case where the
additional neutral Higgs bosons are lighter than $2m_{t}$ as the
phenomenologically interesting region.

We particularly focus on the EW pair production of the additional Higgs
bosons. The cross sections of these processes are simply determined by the
gauge coupling for fixed Higgs boson masses Kanemura and Yuan (2002); Cao  _et
al._ (2004); Belyaev  _et al._ (2008), so that the number of signal events is
determined by decays of the Higgs bosons. This is quite different from single
Higgs productions such as the gluon fusion process and the bottom or top quark
associated process whose cross sections depend on Yukawa couplings. We note
that the single additional Higgs boson production associated with a weak boson
Ellis  _et al._ (1976), e.g., $pp\to H_{2,3}^{0}Z$, is absent at tree level in
the Higgs alignment limit, so that we do not consider these production modes.
Thus, we focus on the EW production, by which we obtain a robust constraint on
the masses of the additional Higgs bosons for given decay branching ratios,
because the production cross section cannot be tuned to be small by taking
model parameters.

###  V.1 EW pair productions of the additional Higgs bosons

We study the pair productions of the additional Higgs bosons. There are the
following 6 subprocesses Eichten  _et al._ (1984); Gunion and Haber (1986);
Djouadi  _et al._ (1999); Cao  _et al._ (2004):

$\displaystyle pp\to H_{2}^{0}H_{3}^{0},~{}~{}pp\to
H_{2}^{0}H^{\pm},~{}~{}pp\to H_{3}^{0}H^{\pm},~{}~{}pp\to H^{+}H^{-}.$ (35)

Figure 3: Cross sections for the EW pair production of the additional Higgs
bosons at LHC with the collision energy of 13 TeV as a function of the
degenerate masses $m_{H_{2}^{0}}(=m_{H_{3}^{0}}=m_{H^{\pm}})$.

Fig. 3 shows the total cross section of the EW pair production at LO in the
unit of pb as a function of the additional Higgs boson mass, where
$m_{H_{2}^{0}}=m_{H_{3}^{0}}=m_{H^{\pm}}$ is assumed. We also show the cross
sections of each subprocess. As we see in this plot, the cross section can be
of ${\cal O}(10-500)$ fb at LHC with 13 TeV for
$m_{H_{2}^{0}},m_{H_{3}^{0}}\lesssim 350$ GeV. Thus, ${\cal O}(10^{3}-10^{5})$
events are already expected even in the current LHC data with the integrated
luminosity of 139 fb-1, and we would expect ${\cal O}(10^{4}-10^{6})$ events
are produced at the HL-LHC with 3 ab-1. 777In this paper, we show the results
for 13 TeV although the ultimate planned center of mass energy at HL-LHC is 14
TeV, where about 10-15 % more events are expected for our signals.

As we discussed in Sec. III, the $\zeta_{f}$ parameters control the branching
ratios of the additional Higgs bosons. Especially, the leptonic branching
ratios are described through the four $R$ parameters; $R$, $R_{\tau}$,
$R^{\pm}$ and $R^{\pm}_{\tau}$, see Eq. (28). The $\zeta_{f}$ dependence of
the model then appears in the signal cross sections including the decay of the
Higgs boson, which can be expressed as

$\displaystyle\sigma(pp\to AB\to i_{A}i_{B})\sim\sigma(pp\to AB)\cdot{\cal
B}(A\to i_{A})\cdot{\cal B}(B\to i_{B}),$ (36)

where $A$ and $B$ are $H_{2}^{0}$, $H_{3}^{0}$ or $H^{\pm}$. Their decay modes
including the subsequent additional Higgs decay modes are indicated by $i_{A}$
and $i_{B}$. For example,
$i_{H_{3}^{0}}=\\{\tau^{+}\tau^{-},b\bar{b},Z^{(\ast)}\tau^{+}\tau^{-},Z^{(\ast)}b\bar{b},W^{\pm(\ast)}\tau^{\mp}\nu,$
$W^{\pm(\ast)}tb,\cdots\\}$ and
$i_{H^{\pm}}=\\{\tau^{\pm}\nu,tb,W^{\pm(\ast)}\tau^{+}\tau^{-},W^{\pm(\ast)}b\bar{b},\cdots\\}$.

When the branching ratios of the additional Higgs bosons into a charged lepton
pair are large, multi-tau signatures from the EW pair production can be used
to probe such models Kanemura  _et al._ (2012); Chun  _et al._ (2015). We will
show that the current multi-lepton searches at LHC provide strong constraints
on the model especially when the additional Higgs bosons dominantly decay into
charged leptons. We will also show the expected constraint on the parameter
space at the HL-LHC by extrapolating these data.

###  V.2 Constraints from multi-lepton searches at LHC

We here consider the constraints on the model by the latest multi-lepton
searches at LHC reported by the ATLAS collaboration Aad  _et al._
(2021).888The corresponding analysis by the CMS collaboration utilizes the
Neural Network approach Tumasyan  _et al._ (2021), which is difficult to
interpret in our model, and thus we do not include it. In this analysis,
events with 4 or more leptons including tau leptons have been searched, which
are classified into 14 signal regions with different combinations of the
requirements on the lepton species, the number of $b$-tagged jets, and the
existence of the lepton pair consistent with a $Z$-boson decay. For these
signal regions, the minimum numbers of isolated electrons, muons,
$\tau$-tagged jets and $b$-tagged jets are required. The following kinematical
cuts are imposed for signal events, i.e., $p_{T}>7$ GeV, $|\eta|<2.47$ for
electrons $p_{T}>5$ GeV, $|\eta|<2.7$ for muons and $p_{T}>20$ GeV,
$|\eta|<2.5$ for jets.

We find that the seven signal regions among them are relevant to constrain our
model, and call the three 3L1T (L: $e$ or $\mu$ and T: $\tau$) signal regions
defined in the ATLAS analysis as SR1-SR3, the three 2L2T regions as SR4-SR6,
and the 5L0T region as SR7. A summary of the definition of the relevant SRs is
shown in Table 2, where the observed and the expected number of SM backgrounds
and that of 95% CL upper bounds for the signal in each SR are also provided.
We have checked that the other seven signal regions provide weaker constraints
than the seven signal regions considered here, and therefore, they are omitted
in the following analysis.

As future prospects at the HL-LHC, we set the constraints on the model
parameters based on the assumption with null contribution from new physics.
The expected upper limit on the number of signal events is then estimated by
using the Poisson distribution, in which we use the expected number of SM
events obtained by extrapolating $N_{\rm exp}^{\rm BG}$ (expected number of SM
events for 139 fb-1) to that given at the HL-LHC, i.e., $N_{\rm exp}^{\rm
BG}\times 3000/139$. For instance for SR6, the upper limit on the signal event
number at 95% CL is given to be about 9.

SR | SR in Ref. Aad  _et al._ (2021) | $N(e,\mu)$ | $N(\tau_{\rm had})$ |  $N$($b$-jets) |  $m_{\rm eff}^{\rm th}$ [GeV] | $N_{\rm obs}^{\rm BG}$ | $N_{\rm exp}^{\rm BG}$ | $S_{\rm obs}^{95}$ | $S_{\rm exp}^{95}$
---|---|---|---|---|---|---|---|---|---
SR1 | ${\tt SR1}^{\rm loose}_{\rm bveto}$ | =3 | $\geq 1$ | $0$ | 600 | 7 | 7.7 | 6.37 | 7.46
SR2 | ${\tt SR1}^{\rm tight}_{\rm bveto}$ | =3 | $\geq 1$ | $0$ | 1000 | 2 | 1.6 | 4.47 | 4.22
SR3 | ${\tt SR1}_{\rm breq}$ | =3 | $\geq 1$ | $\geq 1$ | 1300 | 2 | 2.2 | 4.56 | 4.59
SR4 | ${\tt SR2}^{\rm loose}_{\rm bveto}$ | =2 | $\geq 2$ | $0$ | 600 | 5 | 3.4 | 8.45 | 7.45
SR5 | ${\tt SR2}^{\rm tight}_{\rm bveto}$ | =2 | $\geq 2$ | $0$ | 1000 | 2 | 0.35 | 5.63 | 3.53
SR6 | ${\tt SR2}_{\rm breq}$ | =2 | $\geq 2$ | $\geq 1$ | 1100 | 1 | 0.52 | 4.17 | 3.16
SR7 | ${\tt SR5L}$ | =5 | $\geq 0$ | $\geq 0$ | – | 21 | 12.4 | 17.88 | 9.88
Table 2: Summary of the relevant SRs and corresponding limits for 139 fb-1
quoted from Ref. Aad  _et al._ (2021), where $N_{\rm obs}^{\rm BG}$ ($N_{\rm
exp}^{\rm BG}$) is the observed (expected) number of background events and
$S_{\rm obs}^{95}$ ($S_{\rm exp}^{95}$) is the observed (expected) upper
limits on the number of signal events at 95% CL.

We perform the simulation study by using MadGraph5 Alwall  _et al._ (2011) and
Pythia 8 Sjostrand  _et al._ (2008) with the detector simulation using
Delphes3 de Favereau  _et al._ (2014). We consider the $16$ model points with
different additional Higgs masses labeled with $0\leq n_{2}\leq n_{3}\leq 3$,
corresponding to $180$ GeV $\leq m_{H_{2}^{0}}\leq m_{H_{3}^{0}}\leq 330$ GeV,
as follows:

$\displaystyle m_{H_{2}^{0}}=180+50n_{2}~{}{\rm GeV},\
m_{H_{3}^{0}}=180+50n_{3}~{}{\rm GeV},\
m_{H^{\pm}}=\begin{cases}m_{H_{3}^{0}}&{\rm(heavy}\ H^{\pm}\ {\rm
scenario)}\cr m_{H_{2}^{0}}&{\rm(light}\ H^{\pm}\ {\rm scenario)}\end{cases}.$
(37)

In our simulation, we fix $\zeta_{u}=0.01$, $\zeta_{d}=0.1$ and
$\zeta_{e}=0.5$ for the model points with $|\Delta
m|=m_{H_{3}^{0}}-m_{H_{2}^{0}}\leq 50$ GeV, where we define $\Delta
m=m_{H^{\pm}}-m_{\text{no-deg}}$ with $m_{\text{no-deg}}$ being the mass of
the additional neutral Higgs which is not degenerate with $m_{H^{\pm}}$. For
the model points with $|\Delta m|\geq 100$ GeV, larger $\zeta_{e}$ and
$\zeta_{d}$ values are adopted to achieve $R\sim R_{\tau}\sim 50~{}\%$ to keep
the several decay modes contribute in a similar size.

For different values of $\zeta_{f}$, we estimate the number of events falling
down in a certain SR by multiplying the scaling factors defined as the ratio
of the product of the branching ratios, i.e., ${\cal B}(A\to i_{A}){\cal
B}(B\to i_{B})|_{\zeta_{f}}/{\cal B}(A\to i_{A}){\cal B}(B\to i_{B})|_{\rm
fixed}$ with the numerator (denominator) being the value at arbitrary
$\zeta_{f}$ ($\zeta_{f}$ fixed to be the above values). We note that the
production cross section does not depend on the $\zeta_{f}$ parameters, so
that the scaling factor is needed to be multiplied only to the branching
ratios as explained in the above.

In Fig. 4, we show the region of the parameter space excluded by the LHC data
and flavor experiments at 95 % CL in the $\zeta_{e}$–$\zeta_{d}$ plane for the
case with $(m_{H_{2}^{0}},m_{H_{3}^{0}},m_{H^{\pm}})=(230,280,280)$ GeV (left
panel) and $(230,280,230)$ GeV (right panel). These two mass spectra are the
representative choices for the heavy ${H^{\pm}}$ and the light ${H^{\pm}}$
scenarios. The $\zeta_{u}$ parameter is fixed to be $0.1$ in these plots. The
region below the black-solid (blue-dashed) curves are excluded by using
$S_{\rm obs}^{95}$ ($S_{\rm exp}^{95}$) given in Table 2. The red-dotted curve
shows the expected exclusion of the parameter space at the HL-LHC (3 ab-1). We
note that these excluded regions are obtained by taking into account all the
seven SRs defined in Table 2. Namely, we draw these curves so as to maximize
the area of the excluded region. We find that these constraints are dominantly
determined from SR4, where two isolated leptons and at least two $\tau$-tagged
jets are required with the lowest $m_{\rm eff}$ threshold $m_{\rm eff}^{\rm
th}=600$ GeV. The sensitivity tends to increase as the leptonic (mainly tau)
branching ratio increases, so that it is highly correlated to $R_{\tau}$ which
is the function of $|\zeta_{e}/\zeta_{d}|$ and $|\zeta_{e}/\zeta_{u}|$ as
shown in Fig. 2. Interestingly, even though the light ${H^{\pm}}$ scenario
(shown in the right panel) predicts a larger production cross section than the
heavy $H^{\pm}$ case (shown in the left panel), the constraint from the LHC
data in the former case is weaker. This is because the heavier $H^{\pm}$
provide more $H_{2}^{0}$ via the decay $H^{\pm}\to W^{\pm(*)}H_{2}^{0}$, which
give more tau leptons in the final states. More concretely, we find that the
current LHC data excludes the region with $\zeta_{e}/\zeta_{d}\gtrsim 10$
$(\zeta_{e}/\zeta_{d}\gtrsim 20)$ for the heavy (light) $H^{\pm}$ case as long
as $\zeta_{e}\gtrsim 0.1$ ($\zeta_{e}\gtrsim 3$). These bounds are expected to
be significantly improved at the HL-LHC, where the region with
$\zeta_{e}/\zeta_{d}\gtrsim 2$ can be excluded, shown as the red-dotted
curves.

The constraints from the precision measurements of the flavor observables
discussed in Sec. IV are also overlaid on these panels. The regions excluded
by the measurements of $B\to X_{s}\gamma$, $B_{s}\to\mu\mu$ and the leptonic
$\tau$ decay are depicted in the orange, magenta and cyan regions,
respectively. In addition, those excluded by the current LHC searches for
additional neutral Higgs bosons via the gluon fusion production
($ggH_{2,3}^{0}$) and the bottom quark associated production
$b\bar{b}H_{2,3}^{0}$ with the subsequent $H_{2,3}\to\tau^{+}\tau^{-}$ decay
Aad  _et al._ (2020) are depicted by the black regions. We note that the
flavor constraints are sensitive to the signs (complex phases) of the
$\zeta_{f}$ parameters while the constraints from the LHC multi-lepton
searches are essentially insensitive to those. We see that $B\to X_{s}\gamma$
sets the upper limit on the value of $\zeta_{d}$ to be about 2, which does not
depend on $\zeta_{e}$. On the other hand, the constraint from the single Higgs
productions (black shaded region) excludes the region with larger values of
$\zeta_{e}$ and $\zeta_{d}$, in which the exclusion is almost given only by
the $pp\to b\bar{b}H_{2,3}^{0}\to b\bar{b}\tau^{+}\tau^{-}$ channel. The
region with $\zeta_{e}\gtrsim{\cal O}(100)$ is excluded by the leptonic tau
decay. No further region is excluded by the $B_{s}\to\mu\mu$ data in this
setup. As a summary of this figure, we clarify that except for the bottom-left
region i.e., $\zeta_{e}/\zeta_{d}\lesssim{\cal O}(1)$ with
$\zeta_{d}\lesssim{\cal O}(1)$, most of the parameter region can be excluded
or explored after combining all the constraints (including the expected bound
from the HL-LHC) considered in this section.

Figure 4: Constraint on the parameter space on the $\zeta_{e}$–$\zeta_{d}$
plane in the case with $\zeta_{u}=0.1$ for
$(m_{H_{2}^{0}},m_{H_{3}^{0}},m_{H^{\pm}})=(230,280,280)$ GeV in the left
panel and $(230,280,230)$ GeV in the right panel. The regions below the black-
solid, blue-dashed and red-dotted curves are excluded by the observed LHC Run-
II data, the expected LHC Run-II data and data expected at the HL-LHC,
respectively. The regions shaded by orange, magenta and cyan are excluded by
the constraints from $B\to X_{s}\gamma$, $B_{s}\to\mu\mu$ and the leptonic tau
decay, respectively, while those shaded by black are excluded by the searches
for $gg\to H_{2,3}^{0}\to\tau^{+}\tau^{-}$ and $gg\to b\bar{b}H_{2,3}^{0}\to
b\bar{b}\tau^{+}\tau^{-}$ at the LHC Run-II experiments. The green lines
indicate the case with $\zeta_{d}=\zeta_{e}$ and $\zeta_{d}=\zeta_{e}^{-1}$.

###  V.3 Dependence on the masses of the additional Higgs bosons

We discuss the constraints on the parameter space for the different 16 mass
spectra defined in Eq. (37) in order to see the mass dependence of the
constraints. Figs. 5 and 6 respectively show the results for the heavy
$H^{\pm}$ scenario and the light $H^{\pm}$ scenario. The descriptions for each
curve and shaded region are the same as those explained in Fig. 4.

Figure 5: Constraint on the parameter space on the $\zeta_{e}$–$\zeta_{d}$
plane for $\zeta_{u}=0.1$ in the heavy $H^{\pm}$ scenario
($m_{H^{\pm}}=m_{H_{3}^{0}}\geq m_{H_{2}^{0}}$). The bottom to top panels show
the case with $m_{H_{2}^{0}}=180$, 230, 280 and 330 GeV while the most left
panels show the case with $|\Delta m|(=m_{H_{3}^{0}}-m_{H_{2}^{0}})=0$, and
next panels show the case incrementing $\Delta m$ with a step of 50 GeV (up to
$m_{H_{3}^{0}}=330$ GeV). In each panel, the descriptions of curves and shaded
regions are the same as those given in Fig. 4.

Figure 6: Same as Fig. 5 but for the case with the light $H^{\pm}$ scenario
($m_{H^{\pm}}=m_{H_{2}^{0}}\leq m_{H_{3}^{0}}$).

In Fig. 5, we see that slightly larger regions are excluded for the case with
$|\Delta m|\neq 0$ as compared to that with $|\Delta m|=0$ for
$m_{H_{2}^{0}}=180$ GeV (bottom panels) and 230 GeV (next to bottom panels).
This can be understood by the effect of the $H^{\pm}\to W^{\pm(*)}H_{2}^{0}$
decay with $H_{2}^{0}\to\tau^{+}\tau^{-}$, providing more tau leptons at final
states. We note that the decay of $H^{\pm}\to W^{\pm(*)}H_{2}^{0}$ is
suppressed for larger values of $\zeta_{e}$ and/or $\zeta_{d}$ due to the
enhancement of the $H^{\pm}\to\tau\nu$ and $H^{\pm}\to tb$ modes, so that in
such a region the limit is slightly weaker. For $m_{H_{2}^{0}}=280$ GeV with
$|\Delta m|=50$ GeV (right panel in the second row) and $m_{H_{2}^{0}}=280$
GeV (top panel), no bound is obtained from the current LHC data, which is
simply because of the smaller production cross section as compared with the
case with $|\Delta m|=0$. On the other hand, the HL-LHC sets the limit on
$\zeta_{e}$ and $\zeta_{d}$ for these cases. The behavior of the flavor
constraints is similar to that shown in Fig. 4, while the larger region can be
excluded for the case with smaller $m_{H^{\pm}}$.

In Fig. 6, we see that the smaller area of the parameter region is excluded by
the LHC data as compared with the heavy $H^{\pm}$ scenario shown in Fig. 5.
This is because $H_{3}^{0}$ can decay into not only $Z^{(*)}H_{2}^{0}$ but
also $W^{\pm(*)}H^{\pm}$, where the former produces a tau lepton pair in its
subsequent decay while the latter provides $\tau^{\pm}\nu$. In fact, it is
seen that no bound is obtained for $|\Delta m|\geq 100$ GeV ($|\Delta m|\geq
50$ GeV) for $m_{H_{2}^{0}}=180$ and 230 GeV (280 GeV) from the current LHC
data, while the HL-LHC can exclude a portion of the parameter space in such
cases.

In Appendix A, we show the similar plots for $\zeta_{u}=0$ in Figs. 9 and 10.
The main difference in the collider study from the case with $\zeta_{u}=0.1$
is the absence of $H_{2,3}^{0}\to c\bar{c}$ and the $H^{\pm}\to tb$ decay
being more suppressed. We can see that the above effects slightly make the LHC
bound stronger particularly in the region with smaller $\zeta_{e}$ and
$\zeta_{d}$ values as compared with the case with $\zeta_{u}=0.1$. In
addition, the constraint from $B\to X_{s}\gamma$ disappears in this case, so
that the case with smaller $\zeta_{e}/\zeta_{d}$ values is allowed.

Since the LHC bound is almost described by a line with a constant
$\zeta_{e}/\zeta_{d}$ value, we can extract the value of $\zeta_{e}$ along the
line $\zeta_{e}=\zeta_{d}^{-1}$ as a representative point. In this way, we can
set the upper limit on $\zeta_{e}$ as a function of the mass of the additional
Higgs boson. We here particularly discuss the upper limit on $\zeta_{e}$ in
the Type-X THDM which can be regarded as a special case of the aligned THDM,
i.e., $\zeta_{e}=\zeta_{d}^{-1}=\zeta_{u}^{-1}$, see Eq. (10).

Figure 7: Constraints on the parameter $\zeta_{e}$ in the Type-X THDM
($\zeta_{u}=\zeta_{d}=\zeta_{e}^{-1}$) as a function of the mass of $H^{\pm}$.
The shaded regions are excluded by $B\to X_{s}\gamma$ (blue) and by the
current LHC multi-lepton searches (red) with $\Delta m=0$. The region above
each red (black) curve is excluded by the multi-lepton searches at LHC (the
HL-LHC) with several fixed value of $\Delta m$.

In Fig. 7, we show the constraints on the $\zeta_{e}$ value by the LHC multi-
lepton searches as a function of $m_{H^{\pm}}$ in the Type-X THDM
$(\zeta_{e}=\zeta_{d}^{-1}=\zeta_{u}^{-1})$. The region above the red line is
excluded, where the different lines correspond to the cases with different
$\Delta m$. In particular, we fill out the excluded region by a red color for
$\Delta m=0$. For the same $m_{H^{\pm}}$ value, a larger $\Delta m$ would set
a stronger constraint. To see the difference between the heavy/light $H^{\pm}$
scenarios discussed above, for example, comparing the sensitivity at
$m_{H^{\pm}}=200$ GeV for $\Delta m=50$ GeV and that at $m_{H^{\pm}}=250$ GeV
for $\Delta m=-50$ GeV would be helpful, which share the same $m_{H_{2}^{0}}$
and $m_{H_{3}^{0}}$ values but with different $m_{H^{\pm}}$ values. One can
see the stronger constraint on $\zeta_{e}$ for the larger $m_{H^{\pm}}$ case.
In this figure, the blue shaded region describes the excluded region from the
$B\to X_{s}\gamma$ data. One can see that most of the region is highly
constrained for $\Delta m>0$, and especially most of the parameter space is
already excluded for $\Delta m=150$ GeV when we restrict the masses of the
additional Higgs bosons to be lighter than $2m_{t}$. The future HL-LHC reaches
are shown by the black lines. There are almost no allowed region by combing
the constraints from the HL-LHC and $B\to X_{s}\gamma$. Only exception is the
case with $\Delta m=-100$ GeV, where the $H_{3}^{0}\to H^{\pm}W^{\mp}$ and
$H_{3}^{0}\to H_{2}^{0}Z$ modes dominate for such a small $\zeta_{e}$ case.
The reason why the constraint on $\zeta_{e}$ for $\Delta m=-100$ GeV is weak
compared with the case for $\Delta m=-150$ GeV at $m_{H^{\pm}}=180$ GeV is
that the $H_{3}^{0}\to H_{2}^{0}Z$ mode is relatively small compared with the
$H_{3}^{0}\to H^{\pm}W^{\mp}$ mode by the phase space suppression in the
former case, where the $H_{2}^{0}\to\tau^{+}\tau^{-}$ mode provides an
important contribution to the multi-lepton events. On the other hand, at
$m_{H^{\pm}}=130$ GeV, that for $\Delta m=-100$ GeV is strong compared with
that for $\Delta m=-150$ GeV. This is because the $H_{3}^{0}\to
H^{\pm}W^{\mp}\to(\tau^{+}\nu)(\tau^{-}\nu)$ mode provides an important
contribution instead, where no competition between $\tau\nu$ and $tb$ modes
exists. We conclude that the Type-X THDM can be almost completely probed by
the searches for the additional Higgs bosons at the HL-LHC when their masses
are smaller than the $2m_{t}$. We note again that the constraints from $B\to
X_{s}\gamma$ are sensitive to the complex phases of the $\zeta_{f}$
parameters.

###  V.4 Prospect for mass measurements of additional Higgs bosons

We consider here briefly the prospects of the mass measurements for the
additional Higgs bosons at the HL-LHC. It is important as a basis of more
dedicated study for the property measurements of additional Higgs bosons at
LHC as well as the other future collider experiments including ILC. For
example, the measurements of the CP property of the additional Higgs bosons at
ILC has been discussed assuming the mass measurements are available at LHC
Kanemura  _et al._ (2020). We take the several benchmark points which are
currently not excluded but in future possibly probed at LHC, and demonstrate
the strategy to reconstruct the masses of the additional Higgs bosons. A
summary of the benchmark points taken are shown in Table 3.

scenario |  BPs | $m_{H_{2}^{0}}$ | $m_{H_{3}^{0}}$ | $m_{H^{\pm}}$
---|---|---|---|---
| BP1 | 280 GeV | 180 GeV | 280 GeV
heavy $H^{\pm}$ | BP2 | 280 GeV | 230 GeV | 280 GeV
| BP3 | 230 GeV | 180 GeV | 230 GeV
light $H^{\pm}$ | BP4 | 280 GeV | 180 GeV | 180 GeV
Table 3: Summary of the benchmark points.

We try to access the neutral Higgs boson masses using the
$b\bar{b}\tau^{+}\tau^{-}$ mode, which is expected to be the dominant decay
mode in the available parameter region. It is because the region with too
large ratio of $\zeta_{e}/\zeta_{d}$ is already very constrained by the
current LHC data as shown in the previous sections, which indicates that too
large $R_{\tau}$ is not allowed unless the additional Higgs bosons are
decoupled. As a selection cut we require a event should contain at least two
$b$-tagged jets and at least two tau-tagged jets with $p_{T}>20$ GeV, and
$|\eta|<2.5$. We also require $E\\!\\!\\!/_{T}>50$ GeV. Since a visible
hadronic tau jet carries only a part of the original tau lepton momentum due
to the escaping neutrino momentum, we adopt the collinear approximation
Rainwater  _et al._ (1998) to reconstruct the tau lepton momentum, assuming
the transverse missing momentum $\mathbf{p}\\!\\!\\!/_{T}$ are generated
solely by these neutrino momenta. The explicit procedure is to obtain $c_{1}$
and $c_{2}$ by solving the following relations,

$\displaystyle\mathbf{p}_{\tau_{1}}=(1+c_{1})\mathbf{p}^{\rm vis}_{\tau_{1}},\
\ \ \mathbf{p}_{\tau_{2}}=(1+c_{2})\mathbf{p}^{\rm vis}_{\tau_{2}}\ \ \ \
(c_{1},c_{2}>0),$ (38)
$\displaystyle\mathbf{p}\\!\\!\\!/_{T}=c_{1}\mathbf{p}^{\rm
vis}_{T,\tau_{1}}+c_{2}\mathbf{p}^{\rm vis}_{T,\tau_{2}}.$ (39)

We only accept the events where the above equation has a solution. The larger
the mass of the resonance is, the better this approximation provides the
reconstructed tau momenta, since the momentum carried by the neutrino is
aligned to the visible momentum. With the reconstructed tau momenta and bottom
momenta we can compute $m_{\tau\tau}^{\rm rec}$ and $m_{bb}$ each event. In
Fig. 8, the 2-dimensional $m_{bb}$ vs. $m_{\tau\tau}$ distributions are shown
as a scattering plot in the left panels for BP1 to BP4 from top to bottom. The
dense regions are depicted in red points. In all cases, dense regions are
observed either at the corresponding $(m_{H_{2}^{0}},m_{H_{3}^{0}})$,
$(m_{H_{3}^{0}},m_{H_{2}^{0}})$, and $(m_{H_{2}^{0}},m_{H_{2}^{0}})$ in each
BP.

Figure 8: 2-dimensional $m_{bb}$ vs. $m_{\tau\tau}^{\rm rec}$ distributions
(left panels), $\min(m_{bb},m_{\tau\tau}^{\rm rec})$ distributions (central
panels), and $\max(m_{bb},m_{\tau\tau}^{\rm rec})$ distributions (right
panels) in the $2b2\tau$ events at BP1 - BP4 from top to bottom.

The central and right panels show the corresponding
$\min(m_{bb},m_{\tau\tau}^{\rm rec})$ and $\max(m_{bb},m_{\tau\tau}^{\rm
rec})$ projected distributions, respectively. The expected number of events
for 3 ab-1 at HL-LHC including the effects of the selection cut and the
efficiencies are shown in the plots. The total and the breakdown of the
contributions from the six production modes are plotted. One can see that the
$\min(m_{bb},m_{\tau\tau}^{\rm rec})$ distribution exhibits the peak at the
corresponding $m_{H_{2}^{0}}$ and the $\max(m_{bb},m_{\tau\tau}^{\rm rec})$
distribution does at the $m_{H_{3}^{0}}$ in each BP. For the benchmark points
in the heavy $H^{\pm}$ scenarios (BP1, BP2, BP3), all production modes
contribute to the peak at $m_{H_{2}^{0}}$ in the $\min(m_{bb},m_{\tau\tau})$
distributions. In the $\max(m_{bb},m_{\tau\tau})$ distributions only the
$H_{2}H_{3}$ and $H_{2}H^{\pm}$ production modes contribute to the peak at
$m_{H_{3}^{0}}$ as expected, while all modes contribute to the peak at
$m_{H_{2}^{0}}$. The corresponding $m_{H_{2}^{0}}$ and $m_{H_{3}^{0}}$ values
are indicated by the vertical dashed lines in each plot.

For the light $H^{\pm}$ scenario (BP4), only the $H_{2}H_{3}$ production mode
essentially contributes to these distributions since the decays involving
charged Higgs bosons do not produce events to fulfill the $2b2\tau$ selection
criteria. Therefore, the expected number of the signals is much smaller than
the cases in the heavy $H^{\pm}$ scenario, although the BP4 has the lightest
mass spectrum and the total production cross section is the largest among
those of the four benchmark points. Accessing the charged Higgs mass would be
also possible, for example, following a similar procedure proposed in Ref.
Iguro  _et al._ (2019), although we leave it for a future study.

As we have seen in this section, the masses of the additional Higgs bosons can
be in principle determined using these characteristic peaks in the proposed
distributions in the mass range we consider, whereas it would be rather
difficult when the charged Higgs boson is the lightest among the three
additional Higgs bosons. For quantitative estimation of the precision of the
mass measurements, clearly more dedicated signal and background analyses are
necessary for our benchmark scenarios, as performed in the different occasion
in Ref. Cao  _et al._ (2004). It is beyond the scope of this paper, but we
plan to do it as a future work.

##  VI Conclusions and Dicussions

We have discussed the double-aligned THDM, where the coupling constants of the
discovered Higgs boson to the SM particles are identical to those of the SM
Higgs boson, and the FCNCs do not appear at tree level. This scenario is well
motivated by the scenario of EW baryogenesis, which is compatible with the
current experimental data. In the double-aligned scenario, the additional
Higgs bosons can mainly decay into a fermion pair or a lighter Higgs boson
associated with a gauge boson depending on the $\zeta_{f}$ parameters and the
masses of the additional Higgs bosons. We have discussed the behavior of the
branching ratios in details especially for the CP conserving case. We then
have explicitly shown the critical values of the $\zeta_{f}$ parameters, at
which the sum of the branching ratios to the fermion pairs is 50%. In
addition, we have investigated the current constraints on the parameter space
from the various flavor experiments, i.e., $B\to X_{s}\gamma$,
$B_{s}\to\mu\mu$ and the tau decay, and from the searches in the multi-lepton
final states at the LHC Run-II experiment. For the collider constraints, we
particularly focused on the EW pair production of the additional Higgs bosons
with their masses below $2m_{t}$, whose production cross sections are simply
determined by the masses of the Higgs bosons. It has been found that a large
portion of the parameter space is already excluded by the current LHC data
when the leptonic decays of the additional Higgs bosons are dominant. As a
result, we have found that the remaining parameter space should satisfy
$|\zeta_{d}/\zeta_{e}|\gtrsim{\cal O}(1)$. Our results can be interpreted to
the scenario in the Type-X THDM as a special case. In the case where all the
additional Higgs bosons are degenerate in mass, we have found that the mass
below 190 GeV has been already excluded by the data from the LHC Run-II in the
Type-X THDM. If there is the mass difference the constraint is stronger.
Namely, the charged Higgs mass below 250 GeV (300 GeV) has been excluded for
$\Delta m=100$ GeV (150 GeV). Furthermore, at the HL-LHC we have shown that
most of the parameter region in the Type-X THDM would be explored. Since the
available parameter region tends to predict a significant amount of the
$b\bar{b}$ branching ratios for the additional Higgs bosons, we have
demonstrated the reconstruction of the masses of additional Higgs bosons from
the $b\bar{b}\tau^{+}\tau^{-}$ final state in a few bench mark points.
Extension to the analysis with CP violating phases will be performed
elsewhere.

###### Acknowledgements.

This work was supported, in part, by the Grant-in-Aid on Innovative Areas, the
Ministry of Education, Culture, Sports, Science and Technology, No. 16H06492,
and by the JSPS KAKENHI Grant No. 20H00160 [SK, MT], the Grant-in-Aid for
Scientific Research C, No. 18K03611 [MT] and the Grant-in-Aid for Early-Career
Scientists, No. 19K14714 [KY].

##  Appendix A The results for $\zeta_{u}=0$

In this appendix, we show the plots as Fig. 5 and Fig. 6 in the heavy charged
Higgs scenario and in the light charged Higgs scenario for $\zeta_{u}=0$. Due
to the suppression of the $H_{2,3}\to c\bar{c}$ and $H^{\pm}\to tb$ decay
modes, stronger constraints are in general obtained from the LHC multi-lepton
searches.

Figure 9: The same plots as Fig. 5 but for $\zeta_{u}=0$ in the heavy
$H^{\pm}$ scenario ($m_{H_{2}^{0}}\leq m_{H_{3}^{0}}=m_{H^{\pm}}$).

Figure 10: The same plots as Fig. 6 but for $\zeta_{u}=0$ in the light
$H^{\pm}$ scenario ($m_{H_{2}^{0}}=m_{H^{\pm}}\leq m_{H_{3}^{0}}$).

## References

  * Haber and Kane (1985) H. E. Haber and G. L. Kane, Phys. Rept. 117, 75 (1985).
  * Ma (2001) E. Ma, Phys. Rev. Lett. 86, 2502 (2001), arXiv:hep-ph/0011121 .
  * Zee (1980) A. Zee, Phys. Lett. 93B, 389 (1980), [Erratum: Phys. Lett.95B,461(1980)].
  * Ma (2006) E. Ma, Phys. Rev. D 73, 077301 (2006), arXiv:hep-ph/0601225 .
  * Aoki  _et al._ (2009a) M. Aoki, S. Kanemura, and O. Seto, Phys. Rev. Lett. 102, 051805 (2009a), arXiv:0807.0361 [hep-ph] .
  * Barbieri  _et al._ (2006) R. Barbieri, L. J. Hall, and V. S. Rychkov, Phys. Rev. D 74, 015007 (2006), arXiv:hep-ph/0603188 .
  * Kanemura  _et al._ (2005) S. Kanemura, Y. Okada, and E. Senaha, Phys. Lett. B 606, 361 (2005), arXiv:hep-ph/0411354 .
  * Fromme  _et al._ (2006) L. Fromme, S. J. Huber, and M. Seniuch, JHEP 11, 038 (2006), arXiv:hep-ph/0605242 .
  * Funakubo  _et al._ (1994) K. Funakubo, A. Kakuto, and K. Takenaga, Prog. Theor. Phys. 91, 341 (1994), arXiv:hep-ph/9310267 .
  * Turok and Zadrozny (1991) N. Turok and J. Zadrozny, Nucl. Phys. B 358, 471 (1991).
  * Cline  _et al._ (1996) J. M. Cline, K. Kainulainen, and A. P. Vischer, Phys. Rev. D 54, 2451 (1996), arXiv:hep-ph/9506284 .
  * Modak and Senaha (2019) T. Modak and E. Senaha, Phys. Rev. D 99, 115022 (2019), arXiv:1811.08088 [hep-ph] .
  * Basler  _et al._ (2021) P. Basler, M. Mühlleitner, and J. Müller, (2021), arXiv:2108.03580 [hep-ph] .
  * Glashow and Weinberg (1977) S. L. Glashow and S. Weinberg, Phys. Rev. D 15, 1958 (1977).
  * Grossman (1994) Y. Grossman, Nucl. Phys. B426, 355 (1994), arXiv:hep-ph/9401311 [hep-ph] .
  * Aoki  _et al._ (2009b) M. Aoki, S. Kanemura, K. Tsumura, and K. Yagyu, Phys. Rev. D80, 015017 (2009b), arXiv:0902.4665 [hep-ph] .
  * Barger  _et al._ (1990) V. D. Barger, J. L. Hewett, and R. J. N. Phillips, Phys. Rev. D41, 3421 (1990).
  * Pich and Tuzon (2009) A. Pich and P. Tuzon, Phys. Rev. D 80, 091702 (2009), arXiv:0908.1554 [hep-ph] .
  * Appelquist and Carazzone (1975) T. Appelquist and J. Carazzone, Phys. Rev. D 11, 2856 (1975).
  * Gunion and Haber (2003) J. F. Gunion and H. E. Haber, Phys. Rev. D 67, 075019 (2003), arXiv:hep-ph/0207010 .
  * Davidson and Haber (2005) S. Davidson and H. E. Haber, Phys. Rev. D 72, 035004 (2005), [Erratum: Phys.Rev.D 72, 099902 (2005)], arXiv:hep-ph/0504050 .
  * Aiko  _et al._ (2021) M. Aiko, S. Kanemura, M. Kikuchi, K. Mawatari, K. Sakurai, and K. Yagyu, Nucl. Phys. B 966, 115375 (2021), arXiv:2010.15057 [hep-ph] .
  * Kanemura  _et al._ (2004) S. Kanemura, Y. Okada, E. Senaha, and C. P. Yuan, Phys. Rev. D 70, 115002 (2004), arXiv:hep-ph/0408364 .
  * Kanemura  _et al._ (2015) S. Kanemura, M. Kikuchi, and K. Yagyu, Nucl. Phys. B 896, 80 (2015), arXiv:1502.07716 [hep-ph] .
  * Kanemura and Nagai (2021) S. Kanemura and R. Nagai, (2021), arXiv:2111.12585 [hep-ph] .
  * Braathen and Kanemura (2019) J. Braathen and S. Kanemura, Phys. Lett. B 796, 38 (2019), arXiv:1903.05417 [hep-ph] .
  * Braathen and Kanemura (2020) J. Braathen and S. Kanemura, Eur. Phys. J. C 80, 227 (2020), arXiv:1911.11507 [hep-ph] .
  * Plehn  _et al._ (1996) T. Plehn, M. Spira, and P. M. Zerwas, Nucl. Phys. B 479, 46 (1996), [Erratum: Nucl.Phys.B 531, 655–655 (1998)], arXiv:hep-ph/9603205 .
  * Baur  _et al._ (2002) U. Baur, T. Plehn, and D. L. Rainwater, Phys. Rev. Lett. 89, 151801 (2002), arXiv:hep-ph/0206024 .
  * Baur  _et al._ (2003) U. Baur, T. Plehn, and D. L. Rainwater, Phys. Rev. D 67, 033003 (2003), arXiv:hep-ph/0211224 .
  * Asakawa  _et al._ (2010) E. Asakawa, D. Harada, S. Kanemura, Y. Okada, and K. Tsumura, Phys. Rev. D 82, 115002 (2010), arXiv:1009.4670 [hep-ph] .
  * Barr  _et al._ (2015) A. J. Barr, M. J. Dolan, C. Englert, D. E. Ferreira de Lima, and M. Spannowsky, JHEP 02, 016 (2015), arXiv:1412.7154 [hep-ph] .
  * Papaefstathiou (2015) A. Papaefstathiou, Phys. Rev. D 91, 113016 (2015), arXiv:1504.04621 [hep-ph] .
  * Cao  _et al._ (2017) Q.-H. Cao, G. Li, B. Yan, D.-M. Zhang, and H. Zhang, Phys. Rev. D 96, 095031 (2017), arXiv:1611.09336 [hep-ph] .
  * Gonçalves  _et al._ (2018) D. Gonçalves, T. Han, F. Kling, T. Plehn, and M. Takeuchi, Phys. Rev. D 97, 113004 (2018), arXiv:1802.04319 [hep-ph] .
  * Biekötter  _et al._ (2019) A. Biekötter, D. Gonçalves, T. Plehn, M. Takeuchi, and D. Zerwas, SciPost Phys. 6, 024 (2019), arXiv:1811.08401 [hep-ph] .
  * Abada  _et al._ (2019a) A. Abada _et al._ (FCC), Eur. Phys. J. C 79, 474 (2019a).
  * Abada  _et al._ (2019b) A. Abada _et al._ (FCC), Eur. Phys. J. ST 228, 1109 (2019b).
  * Aihara  _et al._ (2019) H. Aihara _et al._ (ILC), (2019), arXiv:1901.09829 [hep-ex] .
  * Klamka (2021) J. Klamka (CLICdp), in  _European Physical Society Conference on High Energy Physics 2021_ (2021) arXiv:2111.04787 [hep-ex] .
  * Han  _et al._ (2021) T. Han, D. Liu, I. Low, and X. Wang, Phys. Rev. D 103, 013002 (2021), arXiv:2008.12204 [hep-ph] .
  * Chiesa  _et al._ (2020) M. Chiesa, F. Maltoni, L. Mantani, B. Mele, F. Piccinini, and X. Zhao, JHEP 09, 098 (2020), arXiv:2003.13628 [hep-ph] .
  * Kanemura  _et al._ (2020) S. Kanemura, M. Kubota, and K. Yagyu, JHEP 08, 026 (2020), arXiv:2004.03943 [hep-ph] .
  * Andreev  _et al._ (2018) V. Andreev _et al._ (ACME), Nature 562, 355 (2018).
  * Abel  _et al._ (2020) C. Abel _et al._ (nEDM), Phys. Rev. Lett. 124, 081803 (2020), arXiv:2001.11966 [hep-ex] .
  * Cheung  _et al._ (2020) K. Cheung, A. Jueid, Y.-N. Mao, and S. Moretti, Phys. Rev. D 102, 075029 (2020), arXiv:2003.04178 [hep-ph] .
  * Abe  _et al._ (2014) T. Abe, J. Hisano, T. Kitahara, and K. Tobioka, JHEP 01, 106 (2014), [Erratum: JHEP 04, 161 (2016)], arXiv:1311.4704 [hep-ph] .
  * Jung and Pich (2014) M. Jung and A. Pich, JHEP 04, 076 (2014), arXiv:1308.6283 [hep-ph] .
  * Altmannshofer  _et al._ (2020) W. Altmannshofer, S. Gori, N. Hamer, and H. H. Patel, Phys. Rev. D 102, 115042 (2020), arXiv:2009.01258 [hep-ph] .
  * Low  _et al._ (2020) I. Low, N. R. Shah, and X.-P. Wang, (2020), arXiv:2012.00773 [hep-ph] .
  * Enomoto  _et al._ (2021) K. Enomoto, S. Kanemura, and Y. Mura, (2021), arXiv:2111.13079 [hep-ph] .
  * Kanemura  _et al._ (2021) S. Kanemura, M. Kubota, and K. Yagyu, JHEP 04, 144 (2021), arXiv:2101.03702 [hep-ph] .
  * Baer  _et al._ (2013) H. Baer, T. Barklow, K. Fujii, Y. Gao, A. Hoang, S. Kanemura, J. List, H. E. Logan, A. Nomerotski, M. Perelstein, _et al._ , (2013), arXiv:1306.6352 [hep-ph] .
  * Asai  _et al._ (2017) S. Asai, J. Tanaka, Y. Ushiroda, M. Nakao, J. Tian, S. Kanemura, S. Matsumoto, S. Shirai, M. Endo, and M. Kakizaki, (2017), arXiv:1710.08639 [hep-ex] .
  * Fujii  _et al._ (2017) K. Fujii _et al._ , (2017), arXiv:1710.07621 [hep-ex] .
  * Group (2015) C.-S. S. Group, (2015).
  * Bicer  _et al._ (2014) M. Bicer _et al._ (TLEP Design Study Working Group),  _Proceedings, 2013 Community Summer Study on the Future of U.S. Particle Physics: Snowmass on the Mississippi (CSS2013): Minneapolis, MN, USA, July 29-August 6, 2013_ , JHEP 01, 164 (2014), arXiv:1308.6176 [hep-ex] .
  * Kanemura and Yuan (2002) S. Kanemura and C. P. Yuan, Phys. Lett. B 530, 188 (2002), arXiv:hep-ph/0112165 .
  * Cao  _et al._ (2004) Q.-H. Cao, S. Kanemura, and C. P. Yuan, Phys. Rev. D 69, 075008 (2004), arXiv:hep-ph/0311083 .
  * Belyaev  _et al._ (2008) A. Belyaev, Q.-H. Cao, D. Nomura, K. Tobe, and C. P. Yuan, Phys. Rev. Lett. 100, 061801 (2008), arXiv:hep-ph/0609079 .
  * Eberhardt  _et al._ (2021) O. Eberhardt, A. P. n. Martínez, and A. Pich, JHEP 05, 005 (2021), arXiv:2012.09200 [hep-ph] .
  * Barr and Zee (1990) S. M. Barr and A. Zee, Phys. Rev. Lett. 65, 21 (1990), [Erratum: Phys.Rev.Lett. 65, 2920 (1990)].
  * Peskin and Takeuchi (1990) M. E. Peskin and T. Takeuchi, Phys. Rev. Lett. 65, 964 (1990).
  * Bertolini (1986) S. Bertolini, Nucl. Phys. B272, 77 (1986).
  * Peskin and Wells (2001) M. E. Peskin and J. D. Wells, Phys. Rev. D64, 093003 (2001), arXiv:hep-ph/0101342 [hep-ph] .
  * Grimus  _et al._ (2008) W. Grimus, L. Lavoura, O. M. Ogreid, and P. Osland, Nucl. Phys. B801, 81 (2008), arXiv:0802.4353 [hep-ph] .
  * Kanemura  _et al._ (2011) S. Kanemura, Y. Okada, H. Taniguchi, and K. Tsumura, Phys. Lett. B704, 303 (2011), arXiv:1108.3297 [hep-ph] .
  * Aoki  _et al._ (2012) M. Aoki, S. Kanemura, and K. Yagyu, Phys. Rev. D 85, 055007 (2012), arXiv:1110.4625 [hep-ph] .
  * Capdequi Peyranere  _et al._ (1991) M. Capdequi Peyranere, H. E. Haber, and P. Irulegui, Phys. Rev. D 44, 191 (1991).
  * Kanemura (2000) S. Kanemura, Phys. Rev. D 61, 095001 (2000), arXiv:hep-ph/9710237 .
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Alberti  _et al._ (2015) A. Alberti, P. Gambino, K. J. Healey, and S. Nandi, Phys. Rev. Lett. 114, 061802 (2015), arXiv:1411.6560 [hep-ph] .
  * Aoki  _et al._ (2020) S. Aoki _et al._ (Flavour Lattice Averaging Group), Eur. Phys. J. C 80, 113 (2020), arXiv:1902.08191 [hep-lat] .
  * Amhis  _et al._ (2021) Y. S. Amhis _et al._ (HFLAV), Eur. Phys. J. C 81, 226 (2021), arXiv:1909.12524 [hep-ex] .
  * Charles  _et al._ (2015) J. Charles _et al._ , Phys. Rev. D 91, 073007 (2015), arXiv:1501.05013 [hep-ph] .
  * Czakon  _et al._ (2015) M. Czakon, P. Fiedler, T. Huber, M. Misiak, T. Schutzmeier, and M. Steinhauser, JHEP 04, 168 (2015), arXiv:1503.01791 [hep-ph] .
  * Ciuchini  _et al._ (1998) M. Ciuchini, G. Degrassi, P. Gambino, and G. F. Giudice, Nucl. Phys. B 527, 21 (1998), arXiv:hep-ph/9710335 .
  * Misiak and Steinhauser (2017) M. Misiak and M. Steinhauser, Eur. Phys. J. C 77, 201 (2017), arXiv:1702.04571 [hep-ph] .
  * Borzumati and Greub (1998) F. Borzumati and C. Greub, Phys. Rev. D 58, 074004 (1998), arXiv:hep-ph/9802391 .
  * Kagan and Neubert (1999) A. L. Kagan and M. Neubert, Eur. Phys. J. C 7, 5 (1999), arXiv:hep-ph/9805303 .
  * Li  _et al._ (2014) X.-Q. Li, J. Lu, and A. Pich, JHEP 06, 022 (2014), arXiv:1404.5865 [hep-ph] .
  * Bobeth  _et al._ (2014) C. Bobeth, M. Gorbahn, T. Hermann, M. Misiak, E. Stamou, and M. Steinhauser, Phys. Rev. Lett. 112, 101801 (2014), arXiv:1311.0903 [hep-ph] .
  * Krawczyk and Temes (2005) M. Krawczyk and D. Temes, Eur. Phys. J. C 44, 435 (2005), arXiv:hep-ph/0410248 .
  * Ellis  _et al._ (1976) J. R. Ellis, M. K. Gaillard, and D. V. Nanopoulos, Nucl. Phys. B 106, 292 (1976).
  * Eichten  _et al._ (1984) E. Eichten, I. Hinchliffe, K. D. Lane, and C. Quigg, Rev. Mod. Phys. 56, 579 (1984), [Addendum: Rev.Mod.Phys. 58, 1065–1073 (1986)].
  * Gunion and Haber (1986) J. F. Gunion and H. E. Haber, Nucl. Phys. B 278, 449 (1986), [Erratum: Nucl.Phys.B 402, 569–569 (1993)].
  * Djouadi  _et al._ (1999) A. Djouadi, W. Kilian, M. Muhlleitner, and P. M. Zerwas, Eur. Phys. J. C 10, 45 (1999), arXiv:hep-ph/9904287 .
  * Kanemura  _et al._ (2012) S. Kanemura, K. Tsumura, and H. Yokoya, Phys. Rev. D 85, 095001 (2012), arXiv:1111.6089 [hep-ph] .
  * Chun  _et al._ (2015) E. J. Chun, Z. Kang, M. Takeuchi, and Y.-L. S. Tsai, JHEP 11, 099 (2015), arXiv:1507.08067 [hep-ph] .
  * Aad  _et al._ (2021) G. Aad _et al._ (ATLAS),  (2021), 10.1007/JHEP07(2021)167, arXiv:2103.11684 [hep-ex] .
  * Tumasyan  _et al._ (2021) A. Tumasyan _et al._ (CMS), (2021), arXiv:2106.14246 [hep-ex] .
  * Alwall  _et al._ (2011) J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer, and T. Stelzer, JHEP 06, 128 (2011), arXiv:1106.0522 [hep-ph] .
  * Sjostrand  _et al._ (2008) T. Sjostrand, S. Mrenna, and P. Z. Skands, Comput. Phys. Commun. 178, 852 (2008), arXiv:0710.3820 [hep-ph] .
  * de Favereau  _et al._ (2014) J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, and M. Selvaggi (DELPHES 3), JHEP 02, 057 (2014), arXiv:1307.6346 [hep-ex] .
  * Aad  _et al._ (2020) G. Aad _et al._ (ATLAS), Phys. Rev. Lett. 125, 051801 (2020), arXiv:2002.12223 [hep-ex] .
  * Rainwater  _et al._ (1998) D. L. Rainwater, D. Zeppenfeld, and K. Hagiwara, Phys. Rev. D 59, 014037 (1998), arXiv:hep-ph/9808468 .
  * Iguro  _et al._ (2019) S. Iguro, Y. Omura, and M. Takeuchi, JHEP 11, 130 (2019), arXiv:1907.09845 [hep-ph] .