# Flat-band and multi-dimensional fermions in Pb10(PO4)6O4

Botao Fu  College of Physics and Electronic Engineering, Center for
Computational Sciences, Sichuan Normal University, Chengdu, 610068, China  Qin
He  College of Physics and Electronic Engineering, Center for Computational
Sciences, Sichuan Normal University, Chengdu, 610068, China  Xiao-Ping Li
<EMAIL_ADDRESS>School of Physical Science and Technology, Inner Mongolia
University, Hohhot 010021, China

###### Abstract

Employing a combination of first-principles calculations and low-energy
effective models, we present a comprehensive investigation on the electronic
structure of Pb10(PO4)6O4, which exhibits remarkable quasi-one-dimensional
flat-band around the Fermi level that contains novel multi-dimensional
fermions. These flat bands predominantly originate from $p_{x}/p_{y}$ orbital
of the oxygen molecules chain at $4e$ Wyckoff positions, and thus can be well-
captured by a four-band tight-binding model. Furthermore, the abundant crystal
symmetry inherent to Pb10(PO4)6O4 provides an ideal platform for the emergence
of various multi-dimensional fermions, including a 0D four-fold degenerated
Dirac fermion with quadratic dispersion, a 1D quadratic/linear nodal-line
(QNL/LNL) fermion along symmetric $k$-paths, 1D hourglass nodal-line (HNL)
fermion linked to the Dirac fermion, and a 2D symmetry-enforced nodal surface
(NS) found on the $k_{z}$=$\pi$ plane. Moreover, when considering the weak
ferromagnetic order, Pb10(PO4)6O4 transforms into a rare semi-half-metal,
which is characterized by the presence of Dirac fermion and HNL fermion at the
Fermi level for a single spin channel exhibiting 100$\%$ spin polarization.
Our findings reveal the coexistence of flat bands, diverse topological
semimetal states and ferromagnetism within in Pb10(PO4)6O4, which may provide
valuable insights for further exploring intriguing interplay between
superconductivity and exotic electronic states.

##  I Introduction

Recently, attention has turned to a new class of materials with unique
structures, exemplified by the groundbreaking discovery of Cu-doped lead-
apatite (LK-99), Pb10(PO4)6O, which has been reported as the first room-
temperature ambient-pressure superconductorLee  _et al._ (2023a, b).
Unfortunately, as more precise experimental observations accumulateGaristo
(2023), the mystery surrounding superconductivity in LK-99 is gradually
unraveling, and researchers are increasingly inclined to refute the notion of
room-temperature superconductivity in LK-99. Nevertheless, this system is very
interesting in its own electronic structure. The structure is based on the
P63/m-Pb10(PO4)6O structure model as refined by Krivovichev and Burns in
2003Krivovichev and Burns (2003). The ApatitesGriffin (2023) belongs to a
group of materials with the general formula A10(TO4)6X2+x, where A is the
alkaline or rare earth metal, T=Ge, Si, or P, and X=halide, O or OH. The
parent material, Pb10(PO4)6O, can be modulated through various factors, such
as element substitutionsLiu  _et al._ (2011); Baikie  _et al._ (2014); Fleet
_et al._ (2010) and concentration adjustmentsOgawa  _et al._ . Consequently,
the investigation of derivatives of these materials becomes highly
significant.

On the other side, recently, an array of diverse topological quantum states,
including Weyl/Dirac semimetalsXing  _et al._ (2020); Butch  _et al._ (2011);
Zhou  _et al._ (2016); Schimmel  _et al._ (2023); Jia  _et al._ (2019); Li
_et al._ (2021a) and nodal-line semimetalsChen  _et al._ (2017); Muechler  _et
al._ (2019); Li  _et al._ (2021b) have been experimentally observed in several
superconducting systems such as ZrTe5, TaN-series, TaIrTe4, SnSe and YPtBi,
etc. What’s more, the interaction between these emergent topological non-
trivial electronic states and superconductivity has been recognized as a
promising platform for the realization of exotic quasiparticles, holding great
promise for applications in high-speed electronics and topological quantum
computingZhang and Huang (2020); Wu and Wang (2022); Sufyan and Larsson
(2023); Chen (2020).

Previous theoretical work has discussed the possibility of placing O-atom at
any one of the $4e$ Wyckoff Position (WP), which agree that Pb10(PO4)6O is a
insulatorLai  _et al._ (2023); Si and Held (2023). The compound is non-
magnetic and it also has no spin polarization propertiesHao and Fu (2023).
However, strictly speaking, replacing the O atoms are 1/4 occupied with one of
the four equivalent positions, would break the crystal s original high
symmetry, and they have found it plays a crucial role in modulating the
electronic transport properties along the chain direction. T. Kun et al.
predicted that the Pb10(PO4)6O4, with fully occupied $4e$ positions, belongs
to metal and they suggested that the introduction of oxygen vacancies would
lead to the disappearance of the flat bands in the parent materialTao  _et
al._ (2023). While K. Ogawa et al. conducted investigations into the reaction
thermodynamics of Pb10(PO4)6Ox, where $x$ ranges from 0 to 4, thereby
proposing potential synthesis conditions for varying oxygen
concentrationsOgawa  _et al._ . Despite extensive research effortsLai  _et
al._ (2023); Si and Held (2023); Hao and Fu (2023); Tao  _et al._ (2023);
Cabezas-Escares  _et al._ (2023); Liu  _et al._ (2023a); Shen  _et al._
(2023); Liu  _et al._ (2023b), few attention have been payed to the
Pb10(PO4)6O4, especially regarding its underlying topological semimetal state
within the electronic structure and unexplored magnetic properties.

In this work, combining the first-principles calculations and low-energy
effective models, we conduct a comprehensive investigation on the
bandstructure, topological properties and magnetism of Pb10(PO4)6O4 with with
oxygen fully occupied $4e$ positions. First, we reveal origin of quasi-one-
dimensional flat-band around Fermi level by DFT calculation and tight-binding
model. Moreover, via symmetry analysis and effective $k{\cdot}p$ models, we
uncover the coexistence of multi-dimensional fermions around Fermi level,
including 0D quadratic Dirac fermion at A point, 1D QNL/LNL along $\Gamma$A/HK
path, 1D HNL on ALM$\Gamma$ plane, and 2D NS on $k_{z}$=$\pi$ plane. In
addition, the oxygen molecules chains contribute a weak ferromagnetic order,
thus Pb10(PO4)6O4 becomes rare semi-half-metal, simultaneously hosting 100$\%$
spin polarization and quadratic Dirac and HNL fermions at the Fermi level.

Figure 1: Crystal structure of Pb10(PO4)6O4 with bird view (a), top view (b)
and side view (c). (d) The ELF of Pb10(PO4)6O4.

##  II Computational Methods

Structural optimization and electronic structure calculations were carried out
by using Vienna $ab$-$initio$ simulation package (VASP)Kresse and Furthmüller
(1996a, b) with Perdew-Burke-Ernzerhof parameterized generalized gradient
approximation (PBE-GGA) Perdew  _et al._ (1996). The energy cutoff of 500 eV
and a _k_ -point mesh of 5$\times$5$\times$7 are chosen. The ionic relaxations
were performed until the force on each atom was less than 0.01 eV`/`Å and
convergence criterion for the self-consistent electronic minimization loop was
set to 10-6 eV. To analyze the electronic topological properties, the tight-
binding (TB) Hamiltonian was constructed via MagneticTB codeZhang  _et al._
(2022), and surface states were calculated using the iterative Green’s
function method as implemented in the WannierTools packageWu  _et al._
(2018a).

Figure 2: (a)-(c) The band structures of Pb10(PO4)6O4 with different atomic
orbital projections. The 3D view of four bands around Fermi level at $k_{z}=0$
plane in (d) and at $k_{z}=\pi$ plane in (e). (f) The Fermi surface of
Pb10(PO4)6O4. Figure 3: (a) Comparison showing of four-bands of Pb10(PO4)6O4
around Fermi level by DFT and tight-binding (TB) model, where the fitting
parameters of Hamiltonian are chosen as $e_{1}$=-0.224, $t_{1}$=0.165,
$r_{1}$=0.008, $s_{1}$=-0.009 and $s_{2}$=-0.005. The irreducible
representations are marked for each band, where the superscript $\pm$ denotes
the parity eigenvalues. (b) Double degenerated band along $\Gamma$A path. (c)
Double degenerated band along KH path. (d) Hourglass-like band dispersion
along LM path, (e) Hourglass-like band dispersion along LA’ path. The A’ (0,
0, 0.48)located at $\Gamma$A path. (f) 3D view of quadratic Dirac point at A
point on $k_{z}$=$\pi$ plane. (g) Schematically illustrations for the location
of multi-dimensional fermions in the Brillouin zone.

##  III Results and anlysis

###  III.1 Crystal structure and chemical bonding of Pb10(PO4)6O4

The Pb10(PO4)6O4 has the same crystal structure with experimentally
synthesizedLee  _et al._ (2023a) Pb10(PO4)6O with the space group of 176
(P63/m). The optimized lattice constants are $\bm{a}$=$\bm{b}$=10.151 Å and
$\bm{c}$=7.367 Å, slightly larger than those of Pb10(PO4)6O. The schematic
structure is depicted in Fig. 1(a)-(c). In Fig. 1(a), six phosphorus atoms are
situated at the Wyckoff position (WP) of $6h$ (0.5912,0.9643, 0.75) and are
bonded to four nearest neighboring oxygen atoms, forming six [PO4] tetrahedra.
The remaining four oxygen atoms uniformly occupy the $4e$ WP (0, 0, 0.6571),
arranged in a row and neatly positioned on the rhombus of the unit cell. As
shown in Fig. 1(b), the ten lead atoms can be categorized into two groups: the
first group is at the WP of $6h$ (0.2515,0.9882, 0.75), forming a regular
hexagon centered around the oxygen atoms on the vertical edge, the second
group of lead atoms is located at $4f$ (1/3, 2/3, 0.4955), comprising two two
layers as depicted in Fig. 1(c). These lead atoms are surrounded by those
oxygen atoms from the [PO4] tetrahedra, with each lead atom coordinated by six
oxygen atoms, exhibiting bonding characteristics reminiscent of H-MoS2.

The bonding configurations are vividly visualized through calculating the
electron localization function (ELF)Silvi and Savin (1994) in Fig. 1(d). In
the [PO4] tetrahedron, covalent bonds are observed between P and O atoms, with
a bond length of 1.567 Å. Notably, the P atom has no remaining lone pair
electrons, indicating that all its valence electrons are engaged in bonding.
Along the vertical edge of the oxygen atom chain, the intra- and inter-
molecular distances are 1.721 Å and 1.994 Å, respectively. The ELF clearly
reveals the presence of covalent bonds between the O atoms within the
molecule. Additionally, a pair of lone pairs with the characteristics of
$p_{x,y}$-${\pi}^{*}$ anti-bonding orbitals of O2 molecular are clearly
observed. Furthermore, Bader charge analysis indicates that Pb1 and Pb2 atoms
lose 2.59 and 2.60 electrons, respectively, confirming the ionic bonding
between Pb and [PO4], which is consistent with the above analysis.

###  III.2 Quasi-one-dimensional band structure of Pb10(PO4)6O4

The electronic structure, as shown in Fig. 2(a)-(c), reveals that Pb10(PO4)6O4
exhibits metallic properties. Close to the Fermi level, we observe four bands
that are well-separated from the rest bands, hosting pseudo band gaps above
and below them. Interestingly, in the $k_{z}=0$ and $k_{z}=\pi$ planes, these
four bands exhibit strong anisotropic dispersion, where they display nearly
flat dispersion along $k_{x}/k_{y}$ direction, known as “flat bands”, while
showing strong dispersion along the $k_{z}$ direction. To provide a clearer
visualization of these flat bands, we present a 3D view of bands at $k_{z}=0$
and $\pi$ planes in Fig. 2(d) and (e). At $k_{z}=0$ plane, the two conduction
bands (blue color) demonstrate minimal dispersion with a extremely narrow band
width of only 5 meV, indicating a perfect flat band characteristic. Similarly,
the two valence bands (cyan color) also exhibit weak dispersion with a band
width of approximately 60 meV. Moving to the $k_{z}=\pi$ plane, the four bands
become in pairs, intertwined together. Here too, they display very weak
dispersion with a band width of around 23 meV.

To delve deeper into the origins of these distinctive band structures, we
conducted an orbital projection analysis. It’s discovered that the conduction
bands ($\geq$3 eV) mainly arise from the $p$-orbitals of Pb and the
$s$-orbitals of O. Conversely, the valence bands around and below the Fermi
level are primarily contributed by the $s$-orbitals of Pb and the $p$-orbitals
of O. These can be understood as follows: Pb atom loses electron from $6p$
orbitals, resulting in a substantial contribution to the conduction bands.
Simultaneously, the $s$-orbitals of Pb are occupied by electrons, and this
contributes to the formation of the valence bands below the Fermi level (-3 eV
to -1 eV) as depicted in Fig. 2(c). O atom gains electrons, and hence, its
$p$-orbitals mainly contribute to the below and around the Fermi level in Fig.
2(b).

Particularly, through charge density calculations in supplementary materials,
we found that the four energy bands near the Fermi level primarily originate
from the anti-bonding molecular orbital of the oxygen at $4e$ WP positions:

$|1,\pi_{p_{x}}^{*}>,|1,\pi_{p_{y}}^{*}>,|2,\pi_{p_{x}}^{*}>,|2,\pi_{p_{y}}^{*}>,$
(1)

the 1 and 2 refers two O2 molecular, the $\pi_{p_{x,y}}^{*}$ stands for the
anti $\pi$-bond formed by $p_{x,y}$-orbital of O-atoms. Along $z$ direction,
the Hamiltonian is written as,

$\mathcal{H}(k_{z})=e_{0}\tau_{0}\sigma_{0}+t[(1+cos(k_{z})\tau_{x})+sin(k_{z})\tau_{y}]\sigma_{x}.$
(2)

The $\bm{\sigma}$ and $\bm{\tau}$ represent intra- and inter-molecule degrees
of freedom while the $e_{0}$ is the molecular orbital energy level and the $t$
is the coupling strength between two molecules. It gives energy spectrum along
$k_{z}$ as $E_{\pm}=e_{0}\pm 2tcos(k_{z}/2)$ with a gap of $2t$ at $k_{z}$=0.
As $k_{z}$ increases, these bands gradually approach each other, eventually
crossing at $k_{z}$=$\pi$, in consistent with the DFT results in Fig. 3(b).

Previous structural analysis already reveals that these four oxygen atoms form
an effectively one-dimensional oxygen molecule chain along the $z$ direction,
with a large intermolecular distance (9.865 Å), resulting in very weak in-
plane coupling and non-dispersive flat band. As shown in Fig. 2(f), this
quasi-one-dimensional electronic structure is also reflected from the Fermi
surface, where the Fermi surface of the system consists of two parallel
planes, providing perfectly Fermi surface nesting with a vector of (0, 0,
0.517).

###  III.3 Multi-dimensional fermions embedded within the flat band of
Pb10(PO4)6O4

Because electrons near the Fermi level play a dominant role in transport
properties, we will conduct an in-depth analysis of the four energy bands in
close proximity to the Fermi level through through first-principles
calculations, symmetry analysis, and low-energy effective models. A magnified
view of the band structure around the Fermi level is depicted in Fig. 3(a).
Firstly, within the $k_{z}=0$ plane, these four energy bands can be
categorized into two distinct groups. The lower-energy bands represent bonding
states between oxygen molecules, while the higher-energy bands correspond to
anti-bonding states. The separation between these two sets of bands is
dictated by the strength of coupling between molecules along the
$z$-direction, whereas the dispersion within each set is influenced by the in-
plane coupling between O2 molecules. Despite the relatively weak dispersion
observed in the in-plane bands, meticulous calculations unveil that the two
valence bands mutually interact within the energy range of -0.6 eV to -0.5 eV,
resulting in a double degeneracy both at $\Gamma$ and K points. Notably, at
the $\Gamma$ point, there is quadratic dispersion in the $k_{x}$-$k_{y}$
plane, while at K point, a linear dispersion is evident. Upon closer
examination of the two conduction bands, a similar double degeneracy is
observed at both the $\Gamma$ and K points.

Even more intriguingly, these instances of 2-fold degeneracy are not isolated.
Actually, they persist consistently along the $k_{z}$ direction, resulting in
the creation of a quadratic nodal-line (QNL) Yu  _et al._ (2019) along the
$\Gamma$A path and a linear nodal-line (LNL) along the KH path, as illustrated
in Fig. 3(b) and 3(c). As $k_{z}$ evolves, the energy of the conduction bands
decreases while the energy of the valence bands increases until they converge
on the $k_{z}=\pi$ plane. In fact, due to the nonsymmorphic symmetry of 176
space group these bands will maintain a twofold degeneracy throughout the
entire $k_{z}=\pi$ plane, which consequently gives rise to the formation of a
2D nodal surface (NS)Wu  _et al._ (2018b); Khan  _et al._ (2022); Liu  _et
al._ (2021), as shown in Fig. 3(f). Of particular interest is the A point,
where the two doubly nodal surface intersect, forming a fourfold degenerate
Dirac point with quadratic dispersion in the $k_{x}$-$k_{y}$ plane, namely
forming a quadratic Dirac point Wu  _et al._ (2020) (QDP).

Additionally, on the high-symmetric A$\Gamma$ML plane in the Brillouin zone
(BZ), there exists new nodal line structure centered at L point. As
illustrated in Fig. 3(f), along LM and LA’ paths, the 2-fold degenerate
conduction and valence band are lift and inter-exchange partners, leading to a
hourglass-like dispersion relation, creating new double degeneracy point
between middle two-bands. Due to coexistence of space inversion
($\mathcal{P}$) and time reversal symmetry ($\mathcal{T}$) these crossing
points are not isolated, they form a continuous closed curve in the momentum
space. Indeed, it gives rise to a hourglass nodal-line (HNL)Takahashi  _et
al._ (2017); Zheng  _et al._ (2021) centered at L point and linked to A point
as schematically shown in Fig. 3(g). Considering the
$\widetilde{\mathcal{C}}_{6z}$ symmetry, there will be six equivalent HNL
which are linked to the central Dirac point, extending along AL direction,
forming a nodal-network structureNeupane  _et al._ (2016); Fu  _et al._
(2018); Wang  _et al._ (2018) in the extended BZ. As a whole, multi-
dimensional fermions including 0D QDP at A point, 1D QNL/LNL along
$\Gamma$A/KH path, 1D HNL lying on A$\Gamma$ML plane and 2D NS on $k_{z}=\pi$
plane are clearly demonstrated in Pb10(PO4)6O4.

###  III.4  The tight-binding model and low energy effective k$\cdot$p model
of Pb10(PO4)6O4

To unravel the origin of multi-dimensional fermions in Pb10(PO4)6O4, we
construct effective models through symmetry analysis. The lattice of
Pb10(PO4)6O4 belong to a nonsymmorphic space group $P6_{3}/m$ (No. 176) which
hosts various types of symmetry operators such as rotation
$(\mathcal{C}_{3z})$, spatial-inversion ($\mathcal{P}$), mirror
($\mathcal{M}_{z}$) and screw-rotation $(\widetilde{\mathcal{C}}_{2z})$. Based
on the first-principles calculations, we can find that appearance of the four-
fold Dirac point at the $A$ arises from the essential band degeneracy, which
corresponds to the 4D irreducible representations (IRRs) of $A_{2}A_{3}$ [see
Fig. 3(a)]. Under the basis of $A_{2}A_{3}$, the matrix representations of the
generators can be given as

$\displaystyle\widetilde{\mathcal{C}}_{6z}=\frac{\sqrt{3}}{2}\Gamma_{0,3}+\frac{i}{2}\Gamma_{3,3},\thinspace\thinspace\thinspace\mathcal{P}=\Gamma_{0,1},\thinspace\thinspace\thinspace\mathcal{T}=\Gamma_{1,0}\mathcal{K}$
(3)

with $\Gamma_{i,j}=\sigma_{i}\otimes\sigma_{j}$, $\sigma_{0}$ as the $2\times
2$ identity matrix and $\sigma_{i}$ $(i=1,2,3)$ as the Pauli matrix.
$\mathcal{K}$ is a complex conjugate operator. With the standard approach
Zhang  _et al._ (2023), the effective Hamiltonian around the $A$ point
retained to the leading order reads

$\displaystyle\mathcal{H}_{A}$ $\displaystyle=$
$\displaystyle\left[c_{1}+c_{2}(k_{x}^{2}+k_{y}^{2})\right]\Gamma_{0,0}+c_{3}\Gamma_{3,3}k_{z}$
(4) $\displaystyle+$
$\displaystyle\left(\alpha_{1}\Gamma_{+,1}k_{-}^{2}+\alpha_{2}\Gamma_{+,0}k_{+}k_{z}+h.c.\right).$

Here, $\Gamma_{\pm,i}=\sigma_{\pm}\otimes\sigma_{i}$ with
$\sigma_{\pm}=(\sigma_{1}\pm i\sigma_{2})/2$, $k_{\pm}=k_{x}\pm ik_{y}$.
$c_{i=1,2,3}$ is a real parameter, and $\alpha_{i}$ are complex parameters.
Clearly, the obtained Hamiltonian (4) exhibits linear band splitting along the
$k_{z}$ direction and quadratic dispersion in $k_{x}$-$k_{y}$ plane, thus
describes a QDP. Moreover, the QDP hosts a toplogical charge $\mathcal{C}$=0
due to the presence of $\mathcal{P}\mathcal{T}$ symmetry.

Besides the 0D Dirac point, one can observe a doubly degenerate nodal line
along the $\Gamma$A path in Fig. 3(b). The little group of the $\Gamma$A path
is $\mathcal{C}_{6}$, which is an abelian group without any 2D IRRs.
Fortunately, there exists an anti-unitary operator $\mathcal{P}\mathcal{T}$ in
the magnetic co-little group of $\Gamma$A, which can stick two 1D IRRs
together to form a 2D IRRs $\Delta_{3}\Delta_{5}$ or $\Delta_{4}\Delta_{6}$
[see Fig. 3(b)], and thus protects a nodal line existence. We take
$\Delta_{3}\Delta_{5}$ as an example, the matrix representation of symmetry
generators can be expressed as

$\displaystyle\widetilde{\mathcal{C}}_{6z}=-\frac{1}{2}\sigma_{0}+\frac{\sqrt{3}i}{2}\sigma_{3},\thinspace\thinspace\thinspace\thinspace\mathcal{PT}=\sigma_{1}\mathcal{K}.$
(5)

Under the symmetry constraints, the effective Hamiltonian around the $\Gamma$A
path could be derived as

$\displaystyle\mathcal{H}_{{\Gamma}A}$ $\displaystyle=$
$\displaystyle\left[c_{1}+c_{2}k_{z}+c_{3}(k_{x}^{2}+k_{y}^{2})\right]\sigma_{0}$
(6) $\displaystyle+$
$\displaystyle\left(\alpha_{1}\sigma_{+}k_{+}^{2}+h.c.\right),$

where $c_{i}$ are real parameters, and $\alpha_{1}$ is a complex parameters.
The Hamiltonian (6) describes a QNL located at the $\Gamma$A path, which
exhibits linear dispersion along $k_{z}$ and quadratic dispersion normal to
nodal line. For the other 2D IRRs $\Delta_{4}\Delta_{6}$ in the co-little
group of $\Gamma$A, one can check that its effective Hamiltonian shares the
same form as Eq. (6), which also describes a QNL. Besides, applying above
analysis to the double degenerate nodal line along KH path with IRRs
$P_{2}P_{3}$ [see Fig. 3(c)], the corresponding effective Hamiltonian is
written as,

$\displaystyle H_{KH}$ $\displaystyle=$
$\displaystyle\left(c_{1}+c_{2}k_{z}\right)\sigma_{0}+\left(c_{3}k_{x}+c_{4}k_{y}\right)\sigma_{1}$
(7) $\displaystyle+$
$\displaystyle\left(c_{3}k_{y}-c_{4}k_{x}\right)\sigma_{2}$

which describes a LNL as expected.

In addition, there also exists a NS with 2D manifold of reciprocal space due
to the presence of twofold screw-rotational symmetry
$\widetilde{\mathcal{C}}_{2z}=\left\\{\mathcal{C}_{2z}|0,0,\frac{1}{2}\right\\}$
in Pb10(PO4)6O4. It is worth noting that a combination of twofold screw-
rotational symmetry and and time-reversal symmetry $\mathcal{T}$ could
guarantee the appearance of NS on boundary of the Brillouin zone (BZ) Liang
_et al._ (2016). Specifically, one can find that
$(\widetilde{\mathcal{C}}_{2z})^{2}=e^{ik_{z}}$, $\mathcal{T}^{2}=1$ for
spinless case. Therefore, $\widetilde{\mathcal{C}}_{2z}\mathcal{T}$ is anti-
unitary operator and satisfies
$(\widetilde{\mathcal{C}}_{2z}\mathcal{T})^{2}=-1$ on the $k_{z}=\pi$ plane,
which leads to Kramers-like degeneracy of NS. The result is also confirmed by
our DFT calculations.

Since the multi-dimensional fermions including QDP, QNL/LNL, HNL and NS all
originate from the isolated four bands near the Fermi surface, a simple four-
band lattice model is possible to be established, and can be used for further
investigation of the physical properties related to Pb10(PO4)6O4. Based on the
first-principles calculations, we identify the IRRs of four bands near the
Fermi level belong to a single elementary band representations (EBRs)
${}^{1}E^{{}^{\prime\prime}}{}^{2}E^{{}^{\prime\prime}}@2a$ of SG 176 in
topological quantum chemistry Gao  _et al._ (2021); Bradlyn  _et al._ (2017);
Cano  _et al._ (2018). Therefore, we can construct a four-band lattice model
by setting two basis orbitals
$\left(\left|p_{x}\right\rangle,\left|p_{y}\right\rangle\right)$ on $2a$
Wyckoff position in terms of the relationship between the site-symmetry of the
Wyckoff position and the EBRs. The matrix representation of generators of SG
$176$ are given by

$\displaystyle D(\mathcal{C}_{3z})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Gamma_{0,0}-\frac{\sqrt{3}i}{2}\Gamma_{0,2},\thinspace\thinspace\thinspace
D(\widetilde{\mathcal{C}}_{2z})=-\Gamma_{1,0}$ (8) $\displaystyle
D(\mathcal{P})$ $\displaystyle=$
$\displaystyle-\Gamma_{1,0},\thinspace\thinspace\thinspace
D(\mathcal{T})=\Gamma_{0,0}\mathcal{K}.$ (9)

Resultantly, the symmetry allowed lattice model may be written as

$\displaystyle\mathcal{H}_{TB}$ $\displaystyle=$ $\displaystyle
e_{1}+s_{1}\left[\mathrm{cos}\left(k_{x}\right)+\mathrm{cos}\left(k_{y}\right)+\mathrm{cos}\left(k_{x}+k_{y}\right)\right]+2r_{1}\mathrm{cos}\left(k_{z}\right)$
(10) $\displaystyle+$
$\displaystyle\left[\left(\sqrt{3}s_{1}-s_{2}\right)\mathrm{cos}\left(k_{x}\right)+2s_{2}\mathrm{cos}\left(k_{y}\right)-\left(\sqrt{3}s_{1}+s_{2}\right)\mathrm{cos}\left(k_{x}+k_{y}\right)\right]\frac{\Gamma_{0,1}}{2}$
$\displaystyle+$
$\displaystyle\left[2s_{1}\mathrm{cos}\left(k_{y}\right)-\left(\sqrt{3}s_{2}+s_{1}\right)\mathrm{cos}\left(k_{x}\right)+\left(\sqrt{3}s_{2}-s_{1}\right)\mathrm{cos}\left(k_{x}+k_{y}\right)\right]\frac{\Gamma_{0,3}}{2}$
$\displaystyle+$ $\displaystyle
2t_{1}\mathrm{cos}\left(\frac{k_{z}}{2}\right)\Gamma_{1,0}+s_{2}\left[\mathrm{sin}\left(k_{x}\right)+\mathrm{sin}\left(k_{y}\right)-\mathrm{sin}\left(k_{x}+k_{y}\right)\right]\Gamma_{3,2}.$

Here, $e_{1}$ represents the onsite energy, $t_{1}$ and $r_{1}$ correspond to
the nearest neighbor and next neighbor intra-chain hopping along the
$z$-direction, respectively. Additionally, $s_{i}$ denotes the inter-chain
coupling. In Fig. 3(a), we plot the electronic band structure of Eq. (10), one
can observe the established lattice model indeed captures the main feature of
Pb10(PO4)6O4 material with multi-dimensional Fermions.

Figure 4: Spin polarized band structure of Pb10(PO4)6O4. Insets are the
enlarged bands around A and L points.

##  IV semi-half-metal state of Pb10(PO4)6O4 under ferromagnetic ground state

Building on our previous analysis, the low energy electron states of
Pb10(PO4)6O4 are mainly come from the $4e$ O-atoms. These four oxygen atoms
can be viewed as two pairs of interacting oxygen molecules. Given the inherent
paramagnetism of oxygen molecule (2.0 $\mu_{B}$) due to its half-filled 2-fold
$\pi^{*}$ bond, we further explore the underlying magnetic properties of
Pb10(PO4)6O4. The spin polarized calculation is carrier out and we found the
system possesses a total magnetic moment of 2.0 $\mu_{B}$, primarily
originating from O-atom at $4e$. The magnetic moment of O-atom in Pb10(PO4)6O4
is only half of that of isolated oxygen molecule. The reduction in the
magnetic moment could be interpreted from two perspective. On the one hand,
these oxygen atoms have gained electrons (0.4 $e$) from Pb atoms, leading to
electron occupancy exceeding half-filling. On the other hand, there is
electron hopping between two O2 molecules, weaken the electron correlation
effect. Both of these factors contribute to a reduction in the magnetic moment
of the oxygen atoms compared to isolated O2 molecules.

Moreover, by constructing intra-molecular and inter-molecular anti-
ferromagnetic configurations, we discovered that the ferromagnetic
configuration along the $z$ direction is the ground states (with 261 meV lower
than the non-magnetic state and 217 meV lower than the anti-ferromagnetic
state). In Fig. 4, we present the spin-polarized energy bands, revealing a
spin splitting that leads one channel to become a semiconductor while the
other spin channel crosses the Fermi level. This configuration is known as a
half-metal state. What’s particularly intriguing is the presence of completely
spin-polarized flat bands at $k_{z}$=0 plane and topological semimetal states,
including QDP at A point and the HNL around surrounding L point, all located
at the Fermi energy level. This confirms Pb10(PO4)6O4 belongs to the rare
semi-half-metal, which may exhibit immensely strong correlation effects and
topological transport properties.

##  V Conclusions

In summary, our theoretical investigation has unveiled the presence of flat
bands and novel topological semimetal states in the LK-99 derivative. We have
elucidated that the origin of the nearly flat dispersion around the Fermi
level primarily stems from a quasi 1D O-atom chain. Furthermore, the unique
arrangement of O-atoms with nonsymmorphic symmetry plays a crucial role in
generating a 4-fold QDP at the A point, accompanied by the presence of highly
interesting HNL fermions. To comprehensively capture the topological features
of diverse multi-dimensional fermions embedded in Pb10(PO4)6O4, we have
developed both local low-energy k p models and a global four-band tight-
binding model. Our research not only unveils a promising material but also
opens up avenues for exploring the interplay between nontrivial electronic
states, magnetism, and superconductivity within the Pb10(PO4)6O4.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
(NSFC, Grants No. 12304086, No. 12204330), and Dr. B. Fu also the Sichuan
Normal University for financial support (Grant No. 341829001).

## References

  * Lee  _et al._ (2023a) S. Lee, J.-H. Kim, and Y.-W. Kwon, “The first room-temperature ambient-pressure superconductor,” (2023a), arXiv:2307.12008 .
  * Lee  _et al._ (2023b) S. Lee, J. Kim, H.-T. Kim, S. Im, S. An, and K. H. Auh, “Superconductor pb10-xcux(po4)6o showing levitation at room temperature and atmospheric pressure and mechanism,” (2023b), arXiv:2307.12037 .
  * Garisto (2023) D. Garisto, Nature 620, 705 (2023).
  * Krivovichev and Burns (2003) S. V. Krivovichev and P. C. Burns, Z. Krist-Cryst. Mater. 218, 357 (2003).
  * Griffin (2023) S. M. Griffin, “Origin of correlated isolated flat bands in copper-substituted lead phosphate apatite,” (2023), arXiv:2307.16892 .
  * Liu  _et al._ (2011) X. Liu, M. E. Fleet, S. R. Shieh, and Q. He, Phys. Chem. Miner. 38, 397 (2011).
  * Baikie  _et al._ (2014) T. Baikie, M. Schreyer, F. Wei, J. S. Herrin, C. Ferraris, F. Brink, J. Topolska, R. Piltz, J. Price, and T. J. White, Mineral Mag. 78, 325 (2014).
  * Fleet  _et al._ (2010) M. E. Fleet, X. Liu, and S. R. Shieh, Phys. Chem. Miner. 37, 1 (2010).
  * Ogawa  _et al._ (0) K. Ogawa, K. Tolborg, and A. Walsh, ACS Energy Lett. 0, 3941 (0).
  * Xing  _et al._ (2020) Y. Xing, Z. Shao, J. Ge, J. Luo, J. Wang, Z. Zhu, J. Liu, Y. Wang, Z. Zhao, J. Yan, _et al._ , Natl. Sci. Rev. 7, 579 (2020).
  * Butch  _et al._ (2011) N. P. Butch, P. Syers, K. Kirshenbaum, A. P. Hope, and J. Paglione, Phys. Rev. B 84, 220504 (2011).
  * Zhou  _et al._ (2016) Y. Zhou, J. Wu, W. Ning, N. Li, Y. Du, X. Chen, R. Zhang, Z. Chi, X. Wang, X. Zhu, _et al._ , Proc. Natl. Acad. Sci. 113, 2904 (2016).
  * Schimmel  _et al._ (2023) S. Schimmel, Y. Fasano, S. Hoffmann, J. Puig, G. Shipunov, D. Baumann, S. Aswartham, B. Büchner, and C. Hess, arXiv:2302.08968 (2023).
  * Jia  _et al._ (2019) Y.-T. Jia, J.-F. Zhao, S.-J. Zhang, S. Yu, G.-Y. Dai, W.-M. Li, L. Duan, G.-Q. Zhao, X.-C. Wang, X. Zheng, _et al._ , Chin. Phys. Lett. 36, 087401 (2019).
  * Li  _et al._ (2021a) X.-P. Li, K. Deng, B. Fu, Y. Li, D.-S. Ma, J. Han, J. Zhou, S. Zhou, and Y. Yao, Phys. Rev. B 103, L081402 (2021a).
  * Chen  _et al._ (2017) X. Chen, P. Lu, X. Wang, Y. Zhou, C. An, Y. Zhou, C. Xian, H. Gao, Z. Guo, C. Park, B. Hou, K. Peng, X. Zhou, J. Sun, Y. Xiong, Z. Yang, D. Xing, and Y. Zhang, Phys. Rev. B 96, 165123 (2017).
  * Muechler  _et al._ (2019) L. Muechler, Z. Guguchia, J.-C. Orain, J. Nuss, L. M. Schoop, R. Thomale, and F. O. Von Rohr, APL Materials 7 (2019), 10.1063/1.5124242.
  * Li  _et al._ (2021b) X.-P. Li, B. Fu, D.-S. Ma, C. Cui, Z.-M. Yu, and Y. Yao, Phys. Rev. B 103, L161109 (2021b).
  * Zhang and Huang (2020) J. Zhang and G. Q. Huang, J. Phys. Condens. Matter 32, 205702 (2020).
  * Wu and Wang (2022) Z. Wu and Y. Wang, Phys. Rev. B 106, 214510 (2022).
  * Sufyan and Larsson (2023) A. Sufyan and J. A. Larsson, ACS Omega 8, 9607 (2023).
  * Chen (2020) J. Chen, Comp. Mater. Sci. 173, 109466 (2020).
  * Lai  _et al._ (2023) J. Lai, J. Li, P. Liu, Y. Sun, and X.-Q. Chen, J. Mater. Sci. Technol. (2023).
  * Si and Held (2023) L. Si and K. Held, “Electronic structure of the putative room-temperature superconductor pb9cu(po4)6o,” (2023), arXiv:2308.00676 .
  * Hao and Fu (2023) L. Hao and E. Fu, J. Mater. Sci. Technol. (2023).
  * Tao  _et al._ (2023) K. Tao, R. Chen, L. Yang, J. Gao, D. Xue, and C. Jia, “The cu induced ultraflat band in the pb10-xcux(po4)6o4 ($x=0,0.5$),” (2023), arXiv:2308.03218 .
  * Cabezas-Escares  _et al._ (2023) J. Cabezas-Escares, N. Barrera, C. Cardenas, and F. Munoz, arXiv:2308.01135 (2023).
  * Liu  _et al._ (2023a) L. Liu, X. Ren, and J. Hu, arXiv:2308.07345 (2023a).
  * Shen  _et al._ (2023) J. Shen, D. Gaines II, S. Shahabfar, Z. Li, D. Kang, S. Griesemer, A. Salgado-Casanova, T.-c. Liu, C.-T. Chou, Y. Xia, _et al._ , arXiv:2308.07941 (2023).
  * Liu  _et al._ (2023b) R. Liu, T. Guo, J. Lu, J. Ren, and T. Ma, arXiv:2308.08454 (2023b).
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Zhang  _et al._ (2022) Z. Zhang, Z.-M. Yu, G.-B. Liu, and Y. Yao, Comput. Phys. Commun. 270, 108153 (2022).
  * Wu  _et al._ (2018a) Q. Wu, S. Zhang, H.-F. Song, M. Troyer, and A. A. Soluyanov, Comput. Phys. Commun. 224, 405 (2018a).
  * Silvi and Savin (1994) B. Silvi and A. Savin, Nature 371, 683 (1994).
  * Yu  _et al._ (2019) Z.-M. Yu, W. Wu, X.-L. Sheng, Y. X. Zhao, and S. A. Yang, Phys. Rev. B 99, 121106 (2019).
  * Wu  _et al._ (2018b) W. Wu, Y. Liu, S. Li, C. Zhong, Z.-M. Yu, X.-L. Sheng, Y. X. Zhao, and S. A. Yang, Phys. Rev. B 97, 115125 (2018b).
  * Khan  _et al._ (2022) M. R. Khan, K. Bu, J.-T. Wang, and C. Chen, Phys. Rev. B 105, 245152 (2022).
  * Liu  _et al._ (2021) X. Liu, C. Yue, S. V. Erohin, Y. Zhu, A. Joshy, J. Liu, A. M. Sanchez, D. Graf, P. B. Sorokin, Z. Mao, J. Hu, and J. Wei, Nano Lett. 21, 4887 (2021).
  * Wu  _et al._ (2020) W. Wu, Z.-M. Yu, X. Zhou, Y. X. Zhao, and S. A. Yang, Phys. Rev. B 101, 205134 (2020).
  * Takahashi  _et al._ (2017) R. Takahashi, M. Hirayama, and S. Murakami, Phys. Rev. B 96, 155206 (2017).
  * Zheng  _et al._ (2021) B. Zheng, F. Zhan, X. Wu, R. Wang, and J. Fan, Phys. Rev. B 104, L060301 (2021).
  * Neupane  _et al._ (2016) M. Neupane, I. Belopolski, M. M. Hosen, D. S. Sanchez, R. Sankar, M. Szlawska, S.-Y. Xu, K. Dimitri, N. Dhakal, P. Maldonado, P. M. Oppeneer, D. Kaczorowski, F. Chou, M. Z. Hasan, and T. Durakiewicz, Phys. Rev. B 93, 201104 (2016).
  * Fu  _et al._ (2018) B. Fu, X. Fan, D. Ma, C.-C. Liu, and Y. Yao, Phys. Rev. B 98, 075146 (2018).
  * Wang  _et al._ (2018) J.-T. Wang, S. Nie, H. Weng, Y. Kawazoe, and C. Chen, Phys. Rev. Lett. 120, 026402 (2018).
  * Zhang  _et al._ (2023) Z. Zhang, Z.-M. Yu, G.-B. Liu, Z. Li, S. A. Yang, and Y. Yao, Comput. Phys. Commun. 290, 108784 (2023).
  * Liang  _et al._ (2016) Q.-F. Liang, J. Zhou, R. Yu, Z. Wang, and H. Weng, Phys. Rev. B 93, 085427 (2016).
  * Gao  _et al._ (2021) J. Gao, Q. Wu, C. Persson, and Z. Wang, Comput. Phys. Commun. 261, 107760 (2021).
  * Bradlyn  _et al._ (2017) B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Nature 547, 298 (2017).
  * Cano  _et al._ (2018) J. Cano, B. Bradlyn, Z. Wang, L. Elcoro, M. G. Vergniory, C. Felser, M. I. Aroyo, and B. A. Bernevig, Phys. Rev. B 97, 035139 (2018).


aainstitutetext: Max Planck Institute for Physics,
Föhringer Ring 6, 80805 München, Germanybbinstitutetext: Technische
Universität München, Physik-Department,
James-Franck-Strasse 1, 85748 Garching, Germany

# Long-lived particle phenomenology
in the 2HDM+$a$ model

Ulrich Haisch a, b and Luc Schnell<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Higgs decays displaced from the primary interaction vertex represent a
striking experimental signature that is actively searched for by the ATLAS,
CMS and LHCb collaborations. We point out that signals of this type appear in
the context of the 2HDM+$a$ model if the mixing angle $\theta$ of the two CP-
odd weak spin-0 eigenstates is tiny and the dark matter (DM) sector is either
decoupled or kinematically inaccessible. Utilising two suitable benchmark
scenarios, we determine the constraints on the parameter space of the 2HDM+$a$
model that are set by the existing LHC searches for long-lived particles
(LLPs) in Higgs decays. We find that depending on the precise mass spectrum of
the spin-0 states, mixing angles $\theta$ in the ballpark of a few $10^{-8}$
to $10^{-5}$ can be excluded based on LHC Run II data. This finding emphasises
the unique role that searches for displaced signatures can play in
constraining the parameter space of the 2HDM+$a$ model. The ensuing DM
phenomenology is also discussed. In particular, we show that parameter choices
leading to an interesting LLP phenomenology can simultaneously explain the DM
abundance observed in today’s Universe.

††preprint: MPP-2023-27

##  1 Introduction

The two-Higgs-doublet plus pseudoscalar model (2HDM+$a$) Ipek  _et al._
(2014); No (2016); Goncalves  _et al._ (2017); Bauer  _et al._ (2017) has by
now established itself as a pillar of the LHC dark matter (DM) search
programme Aaboud  _et al._ (2018a); Pani and Polesello (2018); Tunney  _et
al._ (2017); Arcadi  _et al._ (2018); Abe  _et al._ (2020); Sirunyan  _et al._
(2019); Aaboud  _et al._ (2019a); Arcadi  _et al._ (2020); Sirunyan  _et al._
(2021a); Aad  _et al._ (2021a, b, c, 2022a); ATL (2021); Robens (2021);
Argyropoulos  _et al._ (2021); Banerjee  _et al._ (2022); Argyropoulos and
Haisch (2022); ATL (2022a); Aad  _et al._ (2022b); CMS (2022a); Arcadi  _et
al._ (2022). It includes a DM candidate in the form of a Dirac fermion which
is a singlet under the Standard Model (SM) gauge group, four 2HDM spin-0
states and an additional CP-odd mediator that is meant to provide the dominant
portal between the dark and the visible sector. Since in models with
pseudoscalar mediators the DM direct detection constraints are weaker compared
to models with scalar mediators, the former models are more attractive from an
astrophysical point of view since they often allow to reproduce the observed
DM relic abundance in a wider parameter space and with less tuning. These
features admit a host of missing transverse momentum ($E_{T}^{\rm miss}$) and
non-$E_{T}^{\rm miss}$ signatures in the 2HDM+$a$ model at colliders which can
and have been consistently compared and combined. See for instance
Argyropoulos  _et al._ (2021); ATL (2022b); Argyropoulos and Haisch (2022) for
such combinations.

Beyond the SM (BSM) scenarios in which the hidden and the visible sectors are
connected through a Higgs portal are also being actively probed for at
colliders. One rather generic feature in such BSM models is the appearance of
new electrically neutral long-lived particles (LLPs) that give rise to
displaced vertex signatures in the LHC detectors — see for example Curtin  _et
al._ (2019); Lee  _et al._ (2019); Alimena  _et al._ (2020) for detailed
reviews of theoretical and experimental aspects of LLPs at the LHC. The main
goal of this article is to point out that besides interesting prompt $E_{T,\rm
miss}$ and non-$E_{T,\rm miss}$ signatures, the 2HDM+$a$ model can also have
an attractive LLP phenomenology. In fact, in the 2HDM+$a$ model the role of
the LLP is played by the additional pseudoscalar $a$, which depending on its
mass can be pair produced efficiently in the decays of both the $125\,{\rm
GeV}$ Higgs and the non-SM CP-even Higgs, i.e. $h\to aa$ and $H\to aa$. To
illustrate the different facets of the LLP phenomenology in the 2HDM+$a$ model
we identify two suitable parameter benchmarks. For these benchmark scenarios
we determine the bounds on the mixing angle $\theta$ of the two CP-odd weak
spin-0 eigenstates as a function of the LLP mass that are set by the existing
LHC searches for displaced Higgs decays Aaij  _et al._ (2017); Aaboud  _et
al._ (2018b, 2019b, 2019c); Aad  _et al._ (2020a, b); Sirunyan  _et al._
(2021b); Tumasyan  _et al._ (2021); Aad  _et al._ (2021d); Tumasyan  _et al._
(2022a, b, c); Aad  _et al._ (2022c, d); Tumasyan  _et al._ (2022d); ATL
(2022c); CMS (2022b). It turns out that depending on the precise mass spectrum
of the spin-0 states, mixing angles $\theta$ from around a few $10^{-8}$ to
about $10^{-5}$ can be excluded with LHC Run II data. To the best of our
knowledge, mixing angles $\theta$ in this range cannot be tested by any other
means, which highlights the special role that LLPs searches play in
constraining the parameter space of the 2HDM+$a$ model. In fact, as we will
further demonstrate, parameter choices that lead to an interesting LLP
phenomenology can in general also correctly predict the measured DM relic
density. The regions of 2HDM+$a$ parameter space singled out in our article
therefore deserve, in our humble opinion, dedicated experimental explorations
in future LHC runs.

This work is structured as follows: in Section 2 we detail the theoretical
ingredients that are relevant in the context of this article. Our general
findings concerning the LLP phenomenology in the 2HDM+$a$ model will be
illustrated in Section 3 by considering two suitable parameter benchmark
scenarios as examples. For these two benchmark choices we derive in Section 4
the constraints that the existing LHC searches for LLPs in displaced Higgs
decays place on the 2HDM+$a$ parameter space. In Section 5 we discuss the
resulting DM phenomenology. Section 6 concludes our work.

##  2 2HDM+$a$ model primer

In order to understand under which circumstances the additional pseudoscalar
$a$ in the 2HDM+$a$ model can be an LLP it is useful to recall its partial
decay modes — further details on the structure of the 2HDM+$a$ model can be
found for instance in Bauer  _et al._ (2017); Abe  _et al._ (2020). In the
alignment limit, i.e. $\cos\left(\beta-\alpha\right)=0$, and choosing for
concreteness the Yukawa sector of the 2HDM+$a$ model to be of type-II, one has
at tree level

$\begin{split}\Gamma\left(a\to\chi\bar{\chi}\right)&=\frac{y_{\chi}^{2}}{8\hskip
0.35565pt\pi}\,m_{a}\,\sqrt{1-\frac{4\hskip
0.35565ptm_{\chi}^{2}}{m_{a}^{2}}}\,\cos^{2}\theta\,,\\\\[5.69054pt]
\Gamma\,(a\to f\bar{f})&=\frac{N_{c}^{f}\eta_{f}^{2}\hskip
1.42262pty_{f}^{2}}{16\hskip 0.35565pt\pi}\,m_{a}\,\sqrt{1-\frac{4\hskip
0.35565ptm_{f}^{2}}{m_{a}^{2}}}\,\sin^{2}\theta\,.\end{split}$ (1)

At the one-loop level the pseudoscalar $a$ can also decay to gauge bosons. The
largest partial decay width is the one to gluon pairs. It takes the form

$\Gamma\left(a\to gg\right)=\frac{\alpha_{s}^{2}}{32\hskip
0.35565pt\pi^{3}\hskip
0.35565ptv^{2}}\,m_{a}^{3}\,\left|\sum_{q=t,b,c}\eta_{q}\,f\left(\frac{4\hskip
0.35565ptm_{q}^{2}}{m_{a}^{2}}\right)\right|^{2}\,\sin^{2}\theta\,,$ (2)

with

$f(z)=z\arctan^{2}\left(\frac{1}{\sqrt{z-1}}\right)\,.$ (3)

Here $m_{a}$ is the mass of the pseudoscalar $a$, $m_{\chi}$ is the mass of
the DM particle, $y_{\chi}$ is the Yukawa coupling of the pseudoscalar $a$ to
a pair of DM particles and $\sin\theta$ quantifies the mixing of the two CP-
odd weak spin-0 eigenstates. Furthermore, $N_{c}^{q}=3$, $N_{c}^{l}=1$,
$\eta_{u}=\cot\beta$, $\eta_{d}=\tan\beta$, $\eta_{l}=\tan\beta$ and
$y_{f}=\sqrt{2}\hskip 0.7113ptm_{f}/v$ with $m_{f}$ the mass of the relevant
SM fermion, $v\simeq 246\,{\rm GeV}$ the Higgs vacuum expectation value (VEV)
and $\alpha_{s}$ the strong coupling constant. From the analytic expressions
(1) and (2) it is evident that the pseudoscalar $a$ can only be long-lived if
$\sin\theta$ is sufficiently small, i.e. $\sin\theta\to 0$, and decays to DM
are strongly suppressed/absent which can be achieved either via decoupling,
i.e. $y_{\chi}\to 0$, or by forbidding the process kinematically, i.e.
$m_{\chi}>m_{a}/2$.

Given the strong suppression of the couplings of the pseudoscalar $a$ to SM
fermions in the limit $\sin\theta\to 0$, the only possibility to produce a
long-lived $a$ is via the decay of heavier spin-0 state $\phi$ into a pair of
such pseudoscalars. In the case that the scalar potential is CP conserving the
$\phi$ has to be a CP-even state which implies that in the 2HDM+$a$ model one
can have both decays of the $125\,{\rm GeV}$ Higgs $h$ and the heavy CP-even
Higgs $H$. The corresponding partial decay widths can be written as

$\Gamma\left(\phi\to aa\right)=\frac{g_{\phi aa}^{2}}{32\hskip
0.35565pt\pi}\,m_{\phi}\,\sqrt{1-\frac{4\hskip
0.35565ptm_{a}^{2}}{m_{\phi}^{2}}}\,,$ (4)

with $\phi=h,H$. For $\sin\theta\simeq 0$ the relevant trilinear couplings are
given by Bauer  _et al._ (2017)

$\begin{split}g_{haa}&\simeq-\frac{2\hskip
0.7113ptv}{m_{h}}\left(\lambda_{P1}\cos^{2}\beta+\lambda_{P2}\sin^{2}\beta\right)\,,\\\\[5.69054pt]
g_{Haa}&\simeq\frac{v}{m_{H}}\sin\left(2\beta\right)\left(\lambda_{P1}-\lambda_{P2}\right)\,,\end{split}$
(5)

where $m_{h}\simeq 125\,{\rm GeV}$ is the mass of the SM-like Higgs, while
$\lambda_{P1}$ and $\lambda_{P2}$ are the quartic couplings that appear in the
2HDM+$a$ scalar potential as follows $P^{2}\hskip
0.7113pt\big{(}\lambda_{P1}\hskip
0.7113ptH_{1}^{\dagger}H_{1}+\lambda_{P2}\hskip
0.7113ptH_{2}^{\dagger}H_{2}\big{)}$ (see for example Bauer  _et al._ (2017);
Abe  _et al._ (2020) for the complete expression of the scalar potential).
Here $P$ denotes the additional pseudoscalar in the weak eigenstate basis
which satisfies $P\simeq a$ for $\sin\theta\simeq 0$.

Figure 1:  Examples of tree-level Feynman diagrams representing $pp\to aa$
production via gluon-gluon-fusion (ggF) Higgs production (left) and $pp\to
l^{+}l^{-}aa$ production in associated $Zh$ production (right) in the 2HDM+$a$
model. The possible decay modes of the pseudoscalar $a$ are not shown. Consult
the text for further details.

The trilinear couplings entering (5) can be constrained phenomenologically. In
the case of $g_{haa}$ one can require that the partial decay width
$\Gamma\left(h\to aa\right)$ does not exceed the total decay width
$\Gamma_{h}$ of the $125\,{\rm GeV}$ Higgs as measured directly at the LHC.
For $m_{a}\ll m_{h}$ this leads to the inequality Argyropoulos and Haisch
(2022)

$\left|g_{haa}\right|\lesssim\sqrt{\frac{32\hskip 0.7113pt\pi\hskip
0.7113pt\Gamma_{h}}{m_{h}}}\simeq 0.94\,,$ (6)

where in the last step we have employed the latest 95% confidence level (CL)
bound of $\Gamma_{h}<1.1\,{\rm GeV}$ that follows from the direct LHC
measurements of the total SM-like Higgs width Sirunyan  _et al._ (2017);
Aaboud  _et al._ (2018c). Inserting the first expression of (5) into (6) then
leads to the following relation

$\left|\lambda_{P1}\cos^{2}\beta+\lambda_{P2}\sin^{2}\beta\right|\lesssim
0.24\,.$ (7)

In the case of $g_{Haa}$ one obtains in a similar fashion

$\left|g_{Haa}\right|\lesssim\sqrt{\frac{32\hskip 0.7113pt\pi\hskip
0.7113pt\Gamma_{H}}{m_{H}}}\simeq 3.2\,,$ (8)

where in the final step we have set the ratio between the total decay width
and the mass of the heavy Higgs to $\Gamma_{H}/m_{H}=10\%$. This choice is
motivated by the observation that for significantly larger ratios different
treatments of the $H$ propagator lead to notable changes in the heavy Higgs
production cross section and thus the LLP signal compared to the case of a
Breit-Wigner propagator with fixed width. Combining the second relation in (5)
with (8) it then follows that

$\left|\frac{v}{m_{H}}\sin\left(2\beta\right)\left(\lambda_{P1}-\lambda_{P2}\right)\right|\lesssim
3.2\,.$ (9)

The above discussion should have shown that in the limit $\sin\theta\to 0$, an
LLP signature can arise in the 2HDM+$a$ model from $h$ or $H$ production
followed by the decay of the intermediate Higgs to a pair of pseudoscalars.
Representative tree-level graphs showing $pp\to aa$ production in ggF Higgs
production (left) and $pp\to l^{+}l^{-}aa$ production in associated Higgs
production (right) that appear in the 2HDM+$a$ model can be found in Figure 1.
Notice that in the former case both the $125\,{\rm GeV}$ Higgs and the heavy
CP-even Higgs contribute in the alignment limit. This is not the case for the
latter process as the $HZZ$ vertex vanishes identically in the limit
$\cos\left(\beta-\alpha\right)\to 0$.

##  3 Parameter benchmarks

Besides the phenomenological bounds (7) and (9) that constrain the sizes of
the couplings $\lambda_{P1}$ and $\lambda_{P2}$, the requirement for the
scalar potential to be bounded from below also restricts the quartic couplings
as well as other parameters of the 2HDM+$a$ model. Assuming that
$\lambda_{P1},\lambda_{P2}>0$ and that $\sin\theta\simeq 0$, one finds two
bounded from below conditions that take the form Abe  _et al._ (2020)

$\lambda_{3}>2\lambda\,,\qquad\lambda_{3}\gtrsim-2\lambda\cot^{2}\left(2\beta\right)\,.$
(10)

Here the parameter $\lambda_{3}$ denotes the usual quartic coupling from the
2HDM scalar potential and $\lambda=m_{h}^{2}/(2\hskip 0.7113ptv^{2})\simeq
0.13$ is the cubic SM Higgs self-coupling. In order to fulfill these relations
and to avoid the tight constraints from Higgs and electroweak precision
physics, we make the following common parameter choices

$\lambda_{3}=0.3\,,\quad\cos\left(\beta-\alpha\right)=0\,,\quad\tan\beta=1\,,\quad
m_{H}=m_{A}=m_{H^{\pm}}\,,\quad y_{\chi}=1\,.$ (11)

We furthermore employ a Yukawa sector of type-II throughout this work.

The first 2HDM+$a$ benchmark scenario that we will study as an example to
illustrate the possible LLP phenomenology in the 2HDM+$a$ model is:

$\big{\\{}\lambda_{P1},\lambda_{P2},m_{\chi}\big{\\}}=\big{\\{}2\cdot
10^{-3},2\cdot 10^{-3},170\,{\rm
GeV}\big{\\}}\,,\quad(\text{benchmark~{}I})\,.$ (12)

We furthermore treat $\sin\theta$ and $m_{a}$ as free parameters but require
that $m_{a}<m_{h}/2$ so that the LLP can be pair produced in the decay of the
$125\,{\rm GeV}$ Higgs boson. The precise value of the common heavy Higgs mass
is irrelevant in such a situation and we simply set it to $m_{H}=600\,{\rm
GeV}$ in benchmark I. Notice that the quartic couplings $\lambda_{P1}$ and
$\lambda_{P2}$ have been chosen such that the constraint (7) is easily
fulfilled. In fact, in the limit $m_{a}\to 0$ the benchmark I parameter
choices imply

$\Gamma_{h}=4.15\,{\rm MeV}\,,$ (13)

a value that is very close to the SM prediction of $\Gamma_{h}^{\rm
SM}=4.07\,{\rm MeV}$ Zyla  _et al._ (2020). The corresponding $h\to aa$
branching ratio is

${\rm BR}\left(h\to aa\right)=1.9\%\,.$ (14)

The proper decay length of the pseudoscalar $a$ for masses in the range
$m_{a}\in[20,60]\,{\rm GeV}$ can be approximated by

$\frac{c\hskip 0.35565pt\tau_{a}}{\rm m}\simeq 4.8\cdot
10^{-12}\,\left(\frac{\rm
GeV}{m_{a}}\right)^{0.9}\,\frac{1}{\sin^{2}\theta}\,,$ (15)

which means that for

$\sin\theta\simeq 4.2\cdot 10^{-7}\,,$ (16)

a pseudoscalar of $m_{a}=40\,{\rm GeV}$ has a proper decay length of around
$1\hskip 0.7113pt{\rm m}$. The result (15) includes higher-order QCD
corrections employing the formulas presented in Appendix A of the paper Haisch
_et al._ (2018) as implemented in Argyropoulos and Haisch (2022b). Also notice
that for the choices (12) and assuming that $\sin\theta\simeq 0$, the
additional 2HDM Higgses are all narrow, i.e. $\Gamma_{H}/m_{H}\simeq 2\%$,
$\Gamma_{A}/m_{A}\simeq 4\%$ and $\Gamma_{H^{\pm}}/m_{H^{\pm}}\simeq 4\%$,
with ${\rm BR}\left(H\to t\bar{t}\right)\simeq 100\%$, ${\rm BR}\left(A\to
t\bar{t}\right)\simeq 100\%$ and ${\rm BR}\left(H^{\pm}\to tb\right)\simeq
100\%$.

In our second 2HDM+$a$ benchmark scenario that leads to an interesting LLP
phenomenology, we consider the following parameter choices

$\big{\\{}\lambda_{P1},\lambda_{P2},m_{\chi}\big{\\}}=\big{\\{}3,0,770\,{\rm
GeV}\big{\\}}\,,\quad(\text{benchmark~{}II})\,.$ (17)

The parameters $\sin\theta$, $m_{H}$ and $m_{a}$ are instead treated as input
with the requirements that $m_{a}>m_{h}/2$ and $m_{a}<m_{H}/2$ so that the LLP
can only be pair produced in the decay $H\to aa$ of the heavy CP-even Higgs
boson $H$. Notice that the values $\lambda_{P1}$ and $\lambda_{P2}$ in (17)
satisfy the constraint (9). Taking for example $m_{H}=600\,{\rm GeV}$ and
$m_{a}=150\,{\rm GeV}$, the total decay width of the heavy CP-even Higgs is
given by

$\Gamma_{H}=22\,{\rm GeV}\,,$ (18)

which implies that $\Gamma_{H}/m_{H}=3.7\%$. The corresponding branching
ratios are

${\rm BR}\left(H\to aa\right)=35\%\,,\qquad{\rm BR}\left(H\to
t\bar{t}\right)=65\%\,,$ (19)

meaning that the decays of the heavy Higgs to two LLPs does not have the
largest branching ratio but that di-top decays are more frequent. Notice that
given the structure of the trilinear coupling $g_{Haa}$ in (5) this feature
will be even more pronounced for heavier CP-even Higgs bosons $H$. In the
range $m_{a}\in[100,300]\,{\rm GeV}$, the proper decay length of the
pseudoscalar $a$ is approximately given by

$\frac{c\hskip 0.35565pt\tau_{a}}{\rm m}\simeq 1.2\cdot
10^{-10}\,\left(\frac{\rm
GeV}{m_{a}}\right)^{1.6}\,\frac{1}{\sin^{2}\theta}\,,$ (20)

where again the results of Haisch  _et al._ (2018); Argyropoulos and Haisch
(2022b) have been used. It follows that for

$\sin\theta\simeq 2.1\cdot 10^{-7}\,,$ (21)

a pseudoscalar of $m_{a}=150\,{\rm GeV}$ has a proper decay length of about
$1\hskip 0.7113pt{\rm m}$. Notice finally that in the case of (12) with
$\sin\theta\simeq 0$, one has $\Gamma_{A}/m_{A}\simeq 4\%$ and
$\Gamma_{H^{\pm}}/m_{H^{\pm}}\simeq 4\%$ with ${\rm BR}\left(A\to
t\bar{t}\right)\simeq 100\%$ and ${\rm BR}\left(H^{\pm}\to tb\right)\simeq
100\%$.

##  4 LLP constraints

At the LHC, searches for displaced Higgs boson decays into LLPs have been
carried out by the ATLAS, CMS and LHCb collaborations in different final
states, covering proper decay lengths from around $10^{-3}\hskip 1.42262pt{\rm
m}$ to $10^{3}\hskip 1.42262pt{\rm m}$ Aaij  _et al._ (2017); Aaboud  _et al._
(2018b, 2019b, 2019c); Aad  _et al._ (2020a, b); Sirunyan  _et al._ (2021b);
Tumasyan  _et al._ (2021); Aad  _et al._ (2021d); Tumasyan  _et al._ (2022a,
b, c); Aad  _et al._ (2022c, d); Tumasyan  _et al._ (2022d); ATL (2022c); CMS
(2022b). The LLP mean decay length determines the search strategies and
reconstruction techniques that are employed — see for instance Section 5 of
the review Argyropoulos  _et al._ (2021) for comprehensive descriptions of the
details of the experimental techniques employed in LHC LLP searches.

We first consider the 2HDM+$a$ benchmark I scenario (12) with $m_{a}<m_{h}/2$.
In this case the pseudoscalar $a$ can be pair produced in the decay of the
$125\,{\rm GeV}$ Higgs boson. Depending on whether the mass of the $a$ is
below or above the bottom-quark threshold, the dominant decay modes of the
pseudoscalar are ${\rm BR}\left(a\to c\bar{c}\right)\simeq 53\%$, ${\rm
BR}\left(a\to\tau^{+}\tau^{-}\right)\simeq 38\%$ and ${\rm BR}\left(a\to
gg\right)\simeq 10\%$ or ${\rm BR}\,(a\to b\bar{b})\simeq 85\%$, ${\rm
BR}\left(a\to c\bar{c}\right)\simeq 4\%$, ${\rm
BR}\left(a\to\tau^{+}\tau^{-}\right)\simeq 7\%$ and ${\rm BR}\left(a\to
gg\right)\simeq 3\%$, respectively. LLP searches that target pseudoscalar pair
production in ggF Higgs or associated $Zh$ production (cf. Figure 1) leading
to multi-jet or four-bottom final states therefore provide the most stringent
constraints. Looking for displaced leptons instead leads to significantly
weaker restrictions because of the small leptonic branching ratios.

In Figure 2 we show an assortment of LLP constraints in the $m_{a}\hskip
0.7113pt$–$\hskip 1.42262pt\sin\theta$ plane that apply in the case of (12).
All limits result from LHC searches that consider ggF Higgs production. The
dotted red exclusion corresponds to the search Aaboud  _et al._ (2019c) that
considers displaced hadronic jets in the ATLAS calorimeter (CM) and the muon
spectrometer (MS) Aaboud  _et al._ (2019b), while the dotted blue constraint
instead results from the ATLAS search Aad  _et al._ (2020a) that utilises the
inner detector (ID) and the MS. These searches use up to $36\,{\rm fb}^{-1}$
and $33\,{\rm fb}^{-1}$ of $\sqrt{s}=13\,{\rm TeV}$ data, respectively. The
dotted green (purple) lines represent an upgrade of the MS (CM) search
strategy to $139\,{\rm fb}^{-1}$ of luminosity collected in LHC Run II. The
corresponding limits are reported in the ATLAS publication Aad  _et al._
(2022c) and Aad  _et al._ (2022d), respectively. The dashed yellow contour is
finally the exclusion that derives from the CMS search Tumasyan  _et al._
(2021) which employs the muon endcap and $137\,{\rm fb}^{-1}$ of
$\sqrt{s}=13\,{\rm TeV}$ data. From the figure it is evident that in the
2HDM+$a$ benchmark I scenario the existing LHC searches for displaced Higgs
decays to hadronic jets allow to exclude values of $\sin\theta$ between around
$10^{-7}$ and $10^{-5}$ with the exact bound depending on the mass of the
pseudoscalar $a$. The excluded parameter space corresponds to proper decay
lengths $c\hskip 0.35565pt\tau_{a}$ in the range from around $59\hskip
1.42262pt{\rm m}$ to $0.08\hskip 0.7113pt{\rm m}$. Notice that given the
smallness of the $h\to aa$ branching ratio $\big{(}$cf. (14)$\big{)}$, our
benchmark I scenario easily evades the present bounds on the undetected or
invisible branching ratios of the $125\,{\rm GeV}$ Higgs ATL (2020) that
amount to 19% and 9%, respectively. In fact, even a possible future high-
luminosity LHC (HL-LHC) upper limit on the invisible branching ratio of the
SM-like Higgs of ${\rm BR}\left(h\to{\rm invisible}\right)<2.5\%$ Cepeda  _et
al._ (2019) would not be stringent enough to test (12) indirectly. This
feature underlines the special role that LLP searches for displaced Higgs
decays can play in testing 2HDM+$a$ models with mixing angles $\theta$ close
to zero.

Figure 2:  95% CL exclusion regions in the $m_{a}\hskip 0.7113pt$–$\hskip
1.42262pt\sin\theta$ plane for the 2HDM+$a$ benchmark I scenario (12). The
dotted red, blue, green and purple lines correspond to the limits following
from the ATLAS searches Aaboud  _et al._ (2019b, c), Aad  _et al._ (2020a),
Aad  _et al._ (2022c) and Aad  _et al._ (2022d), respectively. The dashed
yellow curves instead represent the bound that arises from the CMS search
Tumasyan  _et al._ (2021). The parameter space between the lines is
disfavoured. See main text for further details.

Figure 3:  As Figure 2 but for the two 2HDM+$a$ benchmark II scenarios (17).
The upper (lower) panel depicts the results for $m_{H}=600\,{\rm GeV}$
($m_{H}=1000\,{\rm GeV}$). The dotted red, blue, green and purple lines in the
upper panel correspond to the bounds following from Aaboud  _et al._ (2019b,
c), Aad  _et al._ (2020a), Aad  _et al._ (2022c) and Aad  _et al._ (2022d),
respectively. The dotted red exclusion in the lower panel instead represents
the combination of the ATLAS searches Aaboud  _et al._ (2019b, c); Aad  _et
al._ (2022c, d). The shaded parameter regions are disfavoured. Further
explanations can be found in the main text.

Let us now turn our attention to the benchmark II scenario (17). In this case
the parameters are chosen such that an LLP signal may arise from the prompt
decay of the heavy CP-even Higgs, i.e $H\to aa$, followed by the displaced
decays of the pseudoscalars to a pair of SM fermions $a\to f\bar{f}$ or gluons
$a\to gg$. Given our choice of Yukawa sector and $\tan\beta$, the $a$
dominantly decays to the heaviest SM fermion, which means that depending on
the precise value of its mass either $a\to b\bar{b}$ or $a\to t\bar{t}$
provide the largest rate. To illustrate these two possibilities we consider in
benchmark II the mass combination $m_{H}=600\,{\rm GeV}$ with
$m_{a}\in[50,275]\,{\rm GeV}$ as well as $m_{H}=1000\,{\rm GeV}$ with
$m_{a}\in[50,475]\,{\rm GeV}$. At a centre-of-mass energy of
$\sqrt{s}=13\,{\rm TeV}$ the relevant inclusive heavy Higgs production cross
sections are $\sigma\left(pp\to H\right)\simeq 2.0\,{\rm pb}$ and
$\sigma\left(pp\to H\right)\simeq 0.12\,{\rm pb}$ Hig (2022), respectively.
Notice that in the first case and assuming $m_{a}=150\,{\rm GeV}$, the heavy
Higgs branching ratios (19) as well as ${\rm BR}\,(a\to b\bar{b})\simeq 62\%$,
${\rm BR}\left(a\to c\bar{c}\right)\simeq 3\%$, ${\rm
BR}\left(a\to\tau^{+}\tau^{-}\right)\simeq 7\%$ and ${\rm BR}\left(a\to
gg\right)\simeq 28\%$ apply. In the second case, one has instead ${\rm
BR}\left(H\to aa\right)\simeq 9\%$, ${\rm BR}\left(H\to t\bar{t}\right)\simeq
91\%$ and ${\rm BR}\left(a\to t\bar{t}\right)\simeq 100\%$ for
$m_{a}=400\,{\rm GeV}$.

In the upper panel of Figure 3 we display the relevant 95% CL exclusion
regions in the $m_{a}\hskip 0.7113pt$–$\hskip 1.42262pt\sin\theta$ plane that
apply in the case of the 2HDM+$a$ benchmark II scenario with $m_{H}=600\,{\rm
GeV}$. One observes that taken together the ATLAS searches Aaboud  _et al._
(2019b, c); Aad  _et al._ (2020a, 2022c, 2022d) allow to exclude $\sin\theta$
values between $2\cdot 10^{-8}$ and $2\cdot 10^{-6}$. The corresponding proper
decay lengths $c\hskip 0.35565pt\tau_{a}$ range from $53\hskip 1.42262pt{\rm
m}$ to $0.04\hskip 0.7113pt{\rm m}$. The lower plot in Figure 3 shows the
limits on $\sin\theta$ that a combination of the four ATLAS searches Aaboud
_et al._ (2019b, c); Aad  _et al._ (2022c, d) allow to set in the 2HDM+$a$
benchmark II scenario (17) assuming $m_{H}=1000\,{\rm GeV}$. One observes that
the existing LLP searches can exclude mixing angles for the mass points
$m_{a}=50\,{\rm GeV},150\,{\rm GeV},275\,{\rm GeV}$ and $400\,{\rm GeV}$,
while the 2HDM+$a$ realisation with $m_{a}=475\,{\rm GeV}$ remains untested at
present. For pseudoscalar masses below the top-quark threshold, $\sin\theta$
values between around $4\cdot 10^{-8}$ and $2\cdot 10^{-6}$ are excluded,
whereas for $m_{a}=400\,{\rm GeV}$ mixing parameters in the range of about
$4\cdot 10^{-9}$ and $1\cdot 10^{-9}$ are disfavoured. The excluded parameter
space corresponds to $c\hskip 0.35565pt\tau_{a}$ values ranging from around
$9.7\hskip 0.7113pt{\rm m}$ to $0.06\hskip 0.7113pt{\rm m}$. Notice that the
order of magnitude improvement of the constraint on $\sin\theta$ from the
point $m_{a}=275\,{\rm GeV}$ to $m_{a}=400\,{\rm GeV}$ is readily understood
by recalling that in the former case one has $\Gamma_{a}\simeq 16\,{\rm
MeV}\sin^{2}\theta$, while in the latter case $\Gamma_{a}\simeq 12\,{\rm
GeV}\sin^{2}\theta$ as a result of the open $a\to t\bar{t}$ channel. We add
that improving the limits Aad  _et al._ (2022c, d) by a factor of four would
allow to probe our 2HDM+$a$ benchmark II scenario for $m_{H}=1000\,{\rm GeV}$
and $m_{a}=475\,{\rm GeV}$. A final remark concerns the possibility to search
for the heavy CP-even Higgs in the processes $pp\to H\to t\bar{t}$ or $pp\to
t\bar{t}H\to 4t$. While LHC searches for spin-0 resonances in both di-top
Aaboud  _et al._ (2017, 2018d); Sirunyan  _et al._ (2020a) and four-top
production Sirunyan  _et al._ (2020b); Aad  _et al._ (2020c, 2021e); CMS
(2022c) have been performed, it turns out that the existing searches do not
provide any bound on our 2HDM+$a$ benchmark II model for both $m_{H}=600\,{\rm
GeV}$ and $m_{H}=1000\,{\rm GeV}$ in the $m_{a}$ ranges considered in Figure
3. This finding again stresses the unique role that LLPs searches can play in
the 2HDM+$a$ model in constraining the parameter space.

##  5 Relic density

Figure 4:  Feynman diagrams that lead to DM annihilation via
$\chi\bar{\chi}\to ah$ or $\chi\bar{\chi}\to aH$ (left) and $\chi\bar{\chi}\to
aa$ in the 2HDM+$a$ model. The possible decay modes of the pseudoscalar $a$,
the SM-like Higgs $h$ and the heavy CP-even Higgs $H$ are not shown. Further
details are given in the main text.

In order to understand the physics of standard thermal relic freeze-out in
2HDM+$a$ realisations with $\sin\theta\simeq 0$, we first write the cross
section for annihilation of DM into a final state $X$ as

$\sigma\left(\chi\bar{\chi}\to X\right)v_{\rm
rel}=\sigma^{0}_{X}+\sigma^{1}_{X}\hskip 1.42262ptv_{\rm rel}^{2}\,,$ (22)

where $v_{\rm rel}$ is the relative velocity of the DM pair and the
coefficient $\sigma^{0}_{X}$ ($\sigma^{1}_{X}$) describes the $s\hskip
0.7113pt$-wave ($p\hskip 0.7113pt$-wave) contribution.

In the alignment limit, the possible DM annihilation channels involving a
pseudoscalar $a$ are $\chi\bar{\chi}\to a\to f\bar{f}$, $\chi\bar{\chi}\to
a\to ZH$, $\chi\bar{\chi}\to a\to ah$ and $\chi\bar{\chi}\to a\to aH$ for what
concerns $s$-channel processes and $\chi\bar{\chi}\to aa$ with DM exchange in
the $t$-channel (cf. also Abe  _et al._ (2020)). The annihilation cross
sections (22) of the former two reactions are, however, proportional to
$\sin^{2}\theta$ making them numerically irrelevant in the limit
$\sin\theta\to 0$ unless $m_{a}=m_{\chi}/2$. Such highly tuned solutions to
the DM miracle will not be considered in what follows. Similarly, all DM
annihilation contributions involving the exchange of a heavy pseudoscalar $A$
are suppressed by at least two powers of the sine of the mixing angle
$\theta$, so that only the processes depicted in Figure 4 are relevant for the
calculation of the DM abundance in the context of this work.

The annihilation process $\chi\bar{\chi}\to a\to ah$ proceeds via $s$-wave and
we find for the corresponding coefficient the following analytic result

$\sigma^{0}_{ah}=\frac{y_{\chi}^{2}\hskip
1.42262ptg_{haa}^{2}\cos^{2}\theta}{32\hskip
0.35565pt\pi}\,\sqrt{1-\frac{\left(m_{h}+m_{a}\right)^{2}}{4m_{\chi}^{2}}}\,\sqrt{1-\frac{\left(m_{h}-m_{a}\right)^{2}}{4m_{\chi}^{2}}}\,\frac{v^{2}}{\left(m_{a}^{2}-4m_{\chi}^{2}\right)^{2}+m_{a}^{2}\hskip
0.7113pt\Gamma_{a}^{2}}\,,$ (23)

where the expression for $g_{haa}$ in the limit $\sin\theta\to 0$ can be found
in the first line of (5) and $\Gamma_{a}$ denotes the total decay width of the
pseudoscalar $a$. Since $\sigma^{0}_{ah}\neq 0$ we ignore the $p\hskip
0.7113pt$-wave coefficient $\sigma^{1}_{ah}$ below by setting it to zero. The
result for the $s\hskip 0.7113pt$-wave coefficient $\sigma^{0}_{aH}$
describing DM annihilation through $\chi\bar{\chi}\to a\to aH$ is simply
obtained from (23) by the replacements $g_{haa}\to g_{Haa}$ and $m_{h}\to
m_{H}$.

In the case of $\chi\bar{\chi}\to aa$ the annihilation cross section is
instead $p$-wave suppressed (see Albert  _et al._ (2019) for the calculation
of the $t$-channel contribution in the simplified pseudoscalar DM model) and
the corresponding expansion coefficients take the form $\sigma^{0}_{aa}=0$ and

$\sigma^{1}_{aa}=\frac{y_{\chi}^{4}\hskip 0.7113pt\cos^{4}\theta}{24\hskip
0.35565pt\pi}\,\sqrt{1-\frac{m_{a}^{2}}{m_{\chi}^{2}}}\,\frac{m_{\chi}^{2}\left(m_{a}^{2}-m_{\chi}^{2}\right)^{2}}{\left(m_{a}^{2}-2m_{\chi}^{2}\right)^{4}}\,.$
(24)

Figure 5:  Predicted DM relic abundance in the $m_{a}\hskip 0.35565pt$–$\hskip
1.42262ptm_{\chi}$ plane for the 2HDM+$a$ benchmark I parameter choices (12).
The contour lines indicate the value of $\Omega h^{2}/0.12$, meaning that the
regions below (above) 1 correspond to a DM underabundance (overabundance) in
today’s Universe. For additional details we refer the interested reader to the
main text.

Using the velocity expansion (22) the DM relic density after freeze-out can be
approximated by

$\frac{\Omega h^{2}}{0.12}\simeq\frac{1.6\cdot 10^{-10}\,{\rm
GeV}^{-2}\,x_{f}}{\left\langle\sigma\hskip 0.35565ptv_{\rm
rel}\right\rangle_{f}}\,,\qquad\left\langle\sigma v_{\rm
rel}\right\rangle_{f}=\sum_{X}\left(\sigma_{X}^{0}+\frac{3\hskip
0.7113pt\sigma_{X}^{1}}{x_{f}}\right)\,.$ (25)

Here $x_{f}=m_{\chi}/T_{f}\in[20,30]$ with $T_{f}$ the freeze-out temperature
and the sum over $X$ in principle includes all possible final states. As we
have explained above, for $\sin\theta\simeq 0$ and away from the exceptional
points $m_{a}=m_{\chi}/2$, however, only the channels $X=ah,aH,aa$ are
numerically important. In the limit of heavy DM, i.e. $m_{\chi}\gg
m_{a},m_{h},m_{H}$, the velocity-averaged annihilation cross section at the
freeze-out temperature can be further simplified:

$\left\langle\sigma v_{\rm
rel}\right\rangle_{f}\simeq\frac{y_{\chi}^{2}}{128\hskip 0.7113pt\pi\hskip
0.7113ptm_{\chi}^{2}}\left[\frac{\left(g_{haa}^{2}+g_{Haa}^{2}\right)\hskip
0.35565ptv^{2}}{4\hskip
0.7113ptm_{\chi}^{2}}+\frac{y_{\chi}^{2}}{x_{f}}\right]\,.$ (26)

This approximation shows that the $s$-channel ($t$-channel) contributions to
$\left\langle\sigma v_{\rm rel}\right\rangle_{f}$ scale as $1/m_{\chi}^{4}$
$\left(1/m_{\chi}^{2}\right)$ in the limit of infinitely heavy DM.

Figure 6:  As Figure 5 but for the 2HDM+$a$ benchmark II scenario with
$m_{H}=600\,{\rm GeV}$ (upper panel) and $m_{H}=1000\,{\rm GeV}$ (lower
panel), respectively. Additional details can again be found in the main text.

The formulas (25) and (26) represent useful expressions to estimate $\Omega
h^{2}$. In the case of the benchmark I scenario (12) one has
$g_{haa}^{2}\simeq 6\cdot 10^{-5}$ and $g_{Haa}^{2}=0$, and it is thus a good
approximation to neglect the $s$-channel contributions to $\left\langle\sigma
v_{\rm rel}\right\rangle_{f}$. It follows that

$\frac{\Omega
h^{2}}{0.12}\simeq\frac{0.9}{y_{\chi}^{4}}\left(\frac{x_{f}}{25}\right)^{2}\left(\frac{m_{\chi}}{150\,{\rm
GeV}}\right)^{2}\,.$ (27)

Using $x_{f}\simeq 25$ the relic abundance of $\Omega h^{2}=0.120\pm 0.001$ as
determined by Planck Aghanim  _et al._ (2020) is therefore obtained in the
case of (12) for DM masses $m_{\chi}\simeq 160\,{\rm GeV}$ while for parameter
regions with $m_{\chi}\lesssim 160\,{\rm GeV}$ ($m_{\chi}\gtrsim 160\,{\rm
GeV}$) one expects DM underabundance (overabundance). These expectations agree
quite well with the results of our exact relic calculation that have been
obtained with MadDM Backović  _et al._ (2016) and are shown in Figure 5. In
fact, the exact computation for (12) and $m_{a}=30\,{\rm GeV}$ leads to
$\Omega h^{2}=0.118$, while (27) naively predicts a value that is larger by
around 15%. The observed difference can be traced back to the fact that the
MadDM calculation gives $x_{f}\simeq 21$ in the parameter region of interest
and correctly takes into account the phase-space suppression present in (24)
due to the non-zero values of $m_{a}^{2}/m_{\chi}^{2}$. We add that
$\chi\bar{\chi}\to a\to ah$ annihilation represents a relative contribution to
$\Omega h^{2}$ of less than about 1% in the part of the $m_{a}\hskip
0.35565pt$–$\hskip 1.42262ptm_{\chi}$ plane that is depicted in the figure.
Neglecting the $s$-channel contributions in the approximation (27) is hence
fully justified.

In the case of our 2HDM+$a$ benchmark II parameter scenario (17) the coupling
$g_{haa}$ is no longer small, in fact $g_{haa}^{2}\simeq 35$ and we
furthermore have $g_{Haa}^{2}\neq 0$. On the other hand,
$y_{\chi}^{2}/x_{f}\simeq 0.04$ and thus one can neglect the
$\chi\bar{\chi}\to aa$ contribution to the velocity-averaged annihilation
cross section at the freeze-out temperature (26) to first approximation. Doing
this, one obtains the following simple expression

$\frac{\Omega
h^{2}}{0.12}\simeq\frac{43}{y_{\chi}^{2}\left(g_{haa}^{2}+g_{Haa}^{2}\right)}\,\frac{x_{f}}{25}\left(\frac{m_{\chi}}{800\,{\rm
GeV}}\right)^{4}\,,$ (28)

which approximately describes the resulting DM relic abundance for parameter
choices à la benchmark II. From (28) one hence expects that the correct value
of $\Omega h^{2}$ is realised in the case of (17) for $m_{\chi}\simeq
770\,{\rm GeV}$, while for smaller (larger) DM masses one should have $\Omega
h^{2}\gtrsim 0.12$ $\left(\Omega h^{2}\lesssim 0.12\right)$ if
$m_{H}=600\,{\rm GeV}$ which implies $g_{Haa}^{2}\simeq 1.5$. The results of
the corresponding MadDM computation is displayed in the upper panel of Figure
6. We find that for the parameters (17) together with $m_{a}=150\,{\rm GeV}$
as well as $m_{H}=600\,{\rm GeV}$, the exact calculation predicts $\Omega
h^{2}=0.117$, a value less than 5% below the naive expectation. Numerically,
we furthermore obtain that the relative contribution of $\chi\bar{\chi}\to aa$
to the DM relic density is always below $1\%$ in benchmark II with
$m_{H}=600\,{\rm GeV}$, showing that it can be safely neglected in the
derivation of (28).

The result of our ${\tt MadDM}$ scan in the 2HDM+$a$ benchmark II parameter
scenario with $m_{H}=1000\,{\rm GeV}$ is presented in the lower panel of
Figure 6. It is evident from the plot that in this case the above simplistic
formula is not able to capture the more intricate behaviour of the contours of
constant relic density. In fact, this is not a big surprise because in the
derivation of (28) we have assumed that $m_{\chi}\gg m_{a},m_{h},m_{H}$,
however, one has $m_{\chi}<m_{H}$ in the entire $m_{a}\hskip
0.35565pt$–$\hskip 1.42262ptm_{\chi}$ plane considered. Still the values
$m_{\chi}\in[680,740]\,{\rm GeV}$ of the DM mass that lead to the correct DM
abundance for $m_{a}\in[50,500]\,{\rm GeV}$ are only less than 10% smaller
than what one would expect from (28). Like in the case before, it turns out
that the relative contribution of $s$-channel annihilation amounts to more
than 99% of the predicted values of $\Omega h^{2}$. This shows again that DM
annihilation via $\chi\bar{\chi}\to aa$ is phenomenologically irrelevant in
2HDM+$a$ realisations of the type (17).

Before concluding, we add that DM direct detection experiments do not set
relevant constraints on the 2HDM+$a$ benchmarks (12) and (17) for the mixing
angles $\theta\simeq 0$ necessary to have a long-lived $a$. This is a simple
consequence of the fact that the spin-independent DM-nucleon cross section is
suppressed by both a loop factor and two powers of $\sin\theta$ — see the
recent articles Argyropoulos and Haisch (2022); Arcadi  _et al._ (2018); Bell
_et al._ (2018); Abe  _et al._ (2019); Ertas and Kahlhoefer (2019) for
explicit formulas and further explanations.

##  6 Conclusions

The main lesson that can be learnt from the analytic and numerical results
presented in this work is that LHC searches for displaced Higgs decays can
provide unique constraints on 2HDM+$a$ model realisations. In fact, LLP
signatures appear in the context of the 2HDM+$a$ model if the mixing angle
$\theta$ of the two CP-odd weak spin-0 eigenstates is very small and the DM
sector is either decoupled or kinematically inaccessible. In order to
emphasise this generic finding, we have studied two distinct parameter
benchmarks and explored the sensitivity of the existing LHC LLP searches by
performing parameter scans in the $m_{a}\hskip 0.7113pt$–$\hskip
1.42262pt\sin\theta$ plane. The results of these scans can be found in Figures
2 and 3.

In the benchmark I scenario we have chosen the 2HDM+$a$ parameters such that
the $125\,{\rm GeV}$ Higgs boson gives rise to an LLP signature through its
prompt decay $h\to aa$ followed by the displaced decays of the pseudoscalars
to SM fermions such as $a\to b\bar{b}$ or gluons. One important feature that
is nicely illustrated in our benchmark I scan is that the LLP searches for
displaced hadronic jets that have been performed at LHC Run II can probe
regions of parameter space with mixing angles $\theta$ in the ballpark of
$10^{-7}$ to $10^{-5}$ that are presently not accessible by any other means.
In fact, in our benchmark I scenario the predicted $h\to aa$ branching ratio
turns out to be below the target sensitivity that the HL-LHC is expected to
reach on undetected or invisible decays of the $125\,{\rm GeV}$ Higgs boson.
2HDM+$a$ realisations like (12) are therefore unlikely to be testable
indirectly at the LHC in searches for both prompt $h\to aa\to 4f$ production
or through signatures involving a significant amount of $E_{T,\rm miss}$ such
as $h\to{\rm invisible}$ or mono-jet final states.

The benchmark II parameters were instead chosen such that the LLP signal
arises from the prompt decay $H\to aa$ of the heavy CP-even Higgs followed by
the displaced decays of the pseudoscalars to a pair of SM fermions or gluons.
Depending on the precise value of the LLP mass either the $a\to b\bar{b}$ or
the $a\to t\bar{t}$ mode turns out to provide the largest rate. In order to
illustrate these two distinct possibilities, we have considered in benchmark
II the mass combination $m_{H}=600\,{\rm GeV}$ with $m_{a}\in[50,275]\,{\rm
GeV}$ as well as $m_{H}=1000\,{\rm GeV}$ with $m_{a}\in[50,475]\,{\rm GeV}$.
In the former case we found that the existing LHC searches for displaced heavy
Higgs decays provide stringent constraints on $\theta$, excluding values of
the mixing angle between $2\cdot 10^{-8}$ and $2\cdot 10^{-6}$. The limits on
the benchmark II scenario with $m_{H}=1000\,{\rm GeV}$ turned out to be
noticeably weaker than the bounds for $m_{H}=600\,{\rm GeV}$ due to the order
of magnitude smaller inclusive heavy Higgs production cross section. Still,
for the three mass values $m_{a}=50\,{\rm GeV},150\,{\rm GeV}$ and $275\,{\rm
GeV}$, $\theta$ values between $4\cdot 10^{-8}$ and $2\cdot 10^{-6}$ are
excluded, whereas for $m_{a}=400\,{\rm GeV}$ mixing parameters in the range of
$4\cdot 10^{-9}$ and $1\cdot 10^{-9}$ are disfavoured by the searches for the
displaced heavy Higgs decays performed at LHC Run II. We expect that LLP
searches at LHC Run III and beyond will provide sensitivity to 2HDM+$a$ models
à la (17) with heavy CP-even Higgs masses of the order of $1\,{\rm TeV}$ and
pseudoscalar masses $m_{a}$ above the top-quark threshold. Let us finally note
that in future LHC runs it should also be possible to probe 2HDM+$a$ benchmark
II models with a heavy CP-even Higgs that satisfies $m_{H}>2m_{a}$ and
$m_{H}>2m_{t}$ through spin-0 resonance searches in di-top and four-top
production. A detailed analysis of this issue is however beyond the scope of
this work.

We have furthermore demonstrated that parameter choices that give rise to an
interesting LLP phenomenology can simultaneously explain the observed DM
abundance without excessive tuning. Our relic density scans have been given in
Figures 5 and 6. In our benchmark I scenario, we have seen that due to the
smallness of the $g_{haa}$ coupling the value of $\Omega h^{2}$ is entirely
determined by the $t$-channel process $\chi\bar{\chi}\to aa$ which is $p\hskip
0.7113pt$-wave suppressed. For the range $m_{a}<m_{h}/2$ relevant in benchmark
I, we found that the correct relic abundance can be obtained for DM masses of
around $170\,{\rm GeV}$, rather independently of the precise value of $m_{a}$.
In the case of the benchmark II model instead, only the $s$-channel
annihilation contributions $\chi\bar{\chi}\to a\to ah$ and $\chi\bar{\chi}\to
a\to aH$ that are $s\hskip 0.7113pt$-wave turn out to be phenomenologically
relevant. While the precise value of the DM mass for which $\Omega h^{2}\simeq
0.12$ is achieved depends on the chosen parameters such as $m_{a}$ and
$m_{H}$, we saw that DM masses in the range $m_{\chi}\in[680,770]\,{\rm GeV}$
allow to saturate the measured relic abundance in the part of the parameter
space most relevant for the LLP signals triggered by $H\to aa$.

The results and parameter benchmarks presented in our article should provide a
useful starting point for interpretations of future ATLAS, CMS and LHCb
searches for displaced Higgs decays to hadronic jets in the context of the
2HDM+$a$ model. We therefore encourage and look forward to experimental
explorations in this direction.

###### Acknowledgements.

UH thanks Spyros Argyropoulos for interesting discussions that triggered this
research. His comments on the manuscript and his questions concerning the DM
phenomenology of the model where also very valuable in improving our work. The
research of LS is partially supported by the International Max Planck Research
School (IMPRS) on “Elementary Particle Physics” as well as the Collaborative
Research Center SFB1258.

## References

  * Ipek  _et al._ (2014) S. Ipek, D. McKeen, and A. E. Nelson, Phys. Rev. D 90, 055021 (2014), arXiv:1404.3716 [hep-ph].
  * No (2016) J. M. No, Phys. Rev. D 93, 031701 (2016), arXiv:1509.01110 [hep-ph].
  * Goncalves  _et al._ (2017) D. Goncalves, P. A. N. Machado, and J. M. No, Phys. Rev. D 95, 055027 (2017), arXiv:1611.04593 [hep-ph].
  * Bauer  _et al._ (2017) M. Bauer, U. Haisch, and F. Kahlhoefer, JHEP 05, 138 (2017), arXiv:1701.07427 [hep-ph].
  * Aaboud  _et al._ (2018a) M. Aaboud _et al._ (ATLAS), Eur. Phys. J. C 78, 18 (2018a), arXiv:1710.11412 [hep-ex].
  * Pani and Polesello (2018) P. Pani and G. Polesello, Phys. Dark Univ. 21, 8 (2018), arXiv:1712.03874 [hep-ph].
  * Tunney  _et al._ (2017) P. Tunney, J. M. No, and M. Fairbairn, Phys. Rev. D 96, 095020 (2017), arXiv:1705.09670 [hep-ph].
  * Arcadi  _et al._ (2018) G. Arcadi, M. Lindner, F. S. Queiroz, W. Rodejohann, and S. Vogl, JCAP 03, 042 (2018), arXiv:1711.02110 [hep-ph].
  * Abe  _et al._ (2020) T. Abe _et al._ (LHC Dark Matter Working Group), Phys. Dark Univ. 27, 100351 (2020), arXiv:1810.09420 [hep-ex].
  * Sirunyan  _et al._ (2019) A. M. Sirunyan _et al._ (CMS), Eur. Phys. J. C 79, 280 (2019), arXiv:1811.06562 [hep-ex].
  * Aaboud  _et al._ (2019a) M. Aaboud _et al._ (ATLAS), JHEP 05, 142 (2019a), arXiv:1903.01400 [hep-ex].
  * Arcadi  _et al._ (2020) G. Arcadi, G. Busoni, T. Hugle, and V. T. Tenorth, JHEP 06, 098 (2020), arXiv:2001.10540 [hep-ph].
  * Sirunyan  _et al._ (2021a) A. M. Sirunyan _et al._ (CMS), Eur. Phys. J. C 81, 13 (2021a), [Erratum: Eur. Phys. J. C 81, 333 (2021)], arXiv:2008.04735 [hep-ex].
  * Aad  _et al._ (2021a) G. Aad _et al._ (ATLAS), Eur. Phys. J. C 81, 860 (2021a), arXiv:2011.09308 [hep-ex].
  * Aad  _et al._ (2021b) G. Aad _et al._ (ATLAS), JHEP 10, 013 (2021b), arXiv:2104.13240 [hep-ex].
  * Aad  _et al._ (2021c) G. Aad _et al._ (ATLAS), JHEP 11, 209 (2021c), arXiv:2108.13391 [hep-ex].
  * Aad  _et al._ (2022a) G. Aad _et al._ (ATLAS), Phys. Lett. B 829, 137066 (2022a), arXiv:2111.08372 [hep-ex].
  * ATL (2021) _Combination and summary of ATLAS dark matter searches using $139\,{\rm fb}^{-1}$ of $\sqrt{s}=13\,{\rm TeV}$ $pp$ collision data and interpreted in a two-Higgs-doublet model with a pseudoscalar mediator_, (CERN, Geneva, 2021).
  * Robens (2021) T. Robens, Symmetry 13, 2341 (2021), arXiv:2106.02962 [hep-ph].
  * Argyropoulos  _et al._ (2021) S. Argyropoulos, O. Brandt, and U. Haisch, Symmetry 2021, 13 (2021), arXiv:2109.13597 [hep-ph].
  * Banerjee  _et al._ (2022) S. Banerjee, G. Bélanger, D. Bhatia, B. Fuks, and S. Raychaudhuri, JHEP 07, 111 (2022), arXiv:2110.15391 [hep-ph].
  * Argyropoulos and Haisch (2022) S. Argyropoulos and U. Haisch, SciPost Phys. 13, 007 (2022), arXiv:2202.12631 [hep-ph].
  * ATL (2022a) _Search for dark matter produced in association with a Higgs boson decaying to tau leptons at $\sqrt{s}=13\,{\rm TeV}$ with the ATLAS detector_, (CERN, Geneva, 2022).
  * Aad  _et al._ (2022b) G. Aad _et al._ (ATLAS), (2022b), arXiv:2211.13138 [hep-ex].
  * CMS (2022a) _Search for dark matter in final states with a Higgs boson decaying to a pair of $b$-jets and missing transverse momentum at the HL-LHC_, Tech. Rep. (CERN, Geneva, 2022).
  * Arcadi  _et al._ (2022) G. Arcadi, N. Benincasa, A. Djouadi, and K. Kannike, (2022), arXiv:2212.14788 [hep-ph].
  * ATL (2022b) _Dark matter summary plots for $s$-channel, 2HDM+$a$ and Dark Higgs models_, (CERN, Geneva, 2022).
  * Curtin  _et al._ (2019) D. Curtin _et al._ , Rept. Prog. Phys. 82, 116201 (2019), arXiv:1806.07396 [hep-ph].
  * Lee  _et al._ (2019) L. Lee, C. Ohm, A. Soffer, and T.-T. Yu, Prog. Part. Nucl. Phys. 106, 210 (2019), [Erratum: Prog. Part. Nucl. Phys. 122, 103912 (2022)], arXiv:1810.12602 [hep-ph].
  * Alimena  _et al._ (2020) J. Alimena _et al._ , J. Phys. G 47, 090501 (2020), arXiv:1903.04497 [hep-ex] .
  * Aaij  _et al._ (2017) R. Aaij _et al._ (LHCb), Eur. Phys. J. C 77, 812 (2017), arXiv:1705.07332 [hep-ex].
  * Aaboud  _et al._ (2018b) M. Aaboud _et al._ (ATLAS), JHEP 10, 031 (2018b), arXiv:1806.07355 [hep-ex].
  * Aaboud  _et al._ (2019b) M. Aaboud _et al._ (ATLAS), Phys. Rev. D 99, 052005 (2019b), arXiv:1811.07370 [hep-ex].
  * Aaboud  _et al._ (2019c) M. Aaboud _et al._ (ATLAS), Eur. Phys. J. C 79, 481 (2019c), arXiv:1902.03094 [hep-ex].
  * Aad  _et al._ (2020a) G. Aad _et al._ (ATLAS), Phys. Rev. D 101, 052013 (2020a), arXiv:1911.12575 [hep-ex].
  * Aad  _et al._ (2020b) G. Aad _et al._ (ATLAS), Phys. Rev. D 102, 112006 (2020b), arXiv:2005.12236 [hep-ex].
  * Sirunyan  _et al._ (2021b) A. M. Sirunyan _et al._ (CMS), Phys. Rev. D 104, 012015 (2021b), arXiv:2012.01581 [hep-ex].
  * Tumasyan  _et al._ (2021) A. Tumasyan _et al._ (CMS), Phys. Rev. Lett. 127, 261804 (2021), arXiv:2107.04838 [hep-ex].
  * Aad  _et al._ (2021d) G. Aad _et al._ (ATLAS), JHEP 11, 229 (2021d), arXiv:2107.06092 [hep-ex].
  * Tumasyan  _et al._ (2022a) A. Tumasyan _et al._ (CMS), Eur. Phys. J. C 82, 153 (2022a), arXiv:2110.04809 [hep-ex].
  * Tumasyan  _et al._ (2022b) A. Tumasyan _et al._ (CMS), JHEP 03, 160 (2022b), arXiv:2110.13218 [hep-ex].
  * Tumasyan  _et al._ (2022c) A. Tumasyan _et al._ (CMS), JHEP 04, 062 (2022c), arXiv:2112.13769 [hep-ex].
  * Aad  _et al._ (2022c) G. Aad _et al._ (ATLAS), Phys. Rev. D 106, 032005 (2022c), arXiv:2203.00587 [hep-ex].
  * Aad  _et al._ (2022d) G. Aad _et al._ (ATLAS), JHEP 06, 005 (2022d), arXiv:2203.01009 [hep-ex].
  * Tumasyan  _et al._ (2022d) A. Tumasyan _et al._ (CMS), (2022d), arXiv:2205.08582 [hep-ex].
  * ATL (2022c) _Long-lived particle summary plots for Hidden Sector and Dark Photon models_, (CERN, Geneva, 2022).
  * CMS (2022b) _Long-lived particle summary plots_, (CERN, Geneva, 2022).
  * Sirunyan  _et al._ (2017) A. M. Sirunyan _et al._ (CMS), JHEP 11, 047 (2017), arXiv:1706.09936 [hep-ex].
  * Aaboud  _et al._ (2018c) M. Aaboud _et al._ (ATLAS), Phys. Lett. B 784, 345 (2018c), arXiv:1806.00242 [hep-ex].
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Haisch  _et al._ (2018) U. Haisch, J. F. Kamenik, A. Malinauskas, and M. Spira, JHEP 03, 178 (2018), arXiv:1802.02156 [hep-ph].
  * Argyropoulos and Haisch (2022b) S. Argyropoulos and U. Haisch, 2HDM+$a$ width calculator, (2022b).
  * ATL (2020) _A combination of measurements of Higgs boson production and decay using up to $139\,{\rm fb}^{-1}$ of proton–proton collision data at $\sqrt{s}=13\,{\rm TeV}$ collected with the ATLAS experiment_, (CERN, Geneva, 2020).
  * Cepeda  _et al._ (2019) M. Cepeda _et al._ , CERN Yellow Rep. Monogr. 7, 221 (2019), arXiv:1902.00134 [hep-ph].
  * Hig (2022) _BSM Higgs production cross sections at $\sqrt{s}=13\,{\rm TeV}$ (update in CERN Report4 2016)_, (CERN, Geneva, 2022).
  * Aaboud  _et al._ (2017) M. Aaboud _et al._ (ATLAS), Phys. Rev. Lett. 119, 191803 (2017), arXiv:1707.06025 [hep-ex].
  * Aaboud  _et al._ (2018d) M. Aaboud _et al._ (ATLAS), JHEP 12, 039 (2018d), arXiv:1807.11883 [hep-ex].
  * Sirunyan  _et al._ (2020a) A. M. Sirunyan _et al._ (CMS), JHEP 04, 171 (2020a), [Erratum: JHEP 03, 187 (2022)], arXiv:1908.01115 [hep-ex].
  * Sirunyan  _et al._ (2020b) A. M. Sirunyan _et al._ (CMS), Eur. Phys. J. C 80, 75 (2020b), arXiv:1908.06463 [hep-ex].
  * Aad  _et al._ (2020c) G. Aad _et al._ (ATLAS), Eur. Phys. J. C 80, 1085 (2020c), arXiv:2007.14858 [hep-ex].
  * Aad  _et al._ (2021e) G. Aad _et al._ (ATLAS), JHEP 11, 118 (2021e), arXiv:2106.11683 [hep-ex].
  * CMS (2022c) _Evidence for the simultaneous production of four top quarks in proton-proton collisions at $\sqrt{s}=13\,{\rm TeV}$_, (CERN, Geneva, 2022).
  * Albert  _et al._ (2019) A. Albert _et al._ , Phys. Dark Univ. 26, 100377 (2019), arXiv:1703.05703 [hep-ex].
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron. Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO].
  * Backović  _et al._ (2016) M. Backović, A. Martini, K. Kong, O. Mattelaer, and G. Mohlabeng, AIP Conf. Proc. 1743, 060001 (2016), arXiv:1509.03683 [hep-ph].
  * Bell  _et al._ (2018) N. F. Bell, G. Busoni, and I. W. Sanderson, JCAP 08, 017 (2018), [Erratum: JCAP 01, E01 (2019)], arXiv:1803.01574 [hep-ph].
  * Abe  _et al._ (2019) T. Abe, M. Fujiwara, and J. Hisano, JHEP 02, 028 (2019), arXiv:1810.01039 [hep-ph].
  * Ertas and Kahlhoefer (2019) F. Ertas and F. Kahlhoefer, JHEP 06, 052 (2019), arXiv:1902.11070 [hep-ph].