Semi-dark Higgs decays: sweeping the Higgs neutrino floor

We study exotic Higgs decays $h \to Z X$, with $X$ an invisible beyond the Standard Model (SM) particle, resulting in a semi-dark final state. Such exotic Higgs decays may occur in theories of axion-like-particles (ALPs), dark photons or pseudoscalar mediators between the SM and dark matter. The SM process $h\to Z\nu\bar{\nu}$ represents an irreducible"neutrino floor"background to these new physics searches, providing also a target experimental sensitivity for them. We analyze $h \to Z + \text{invisible}$ searches at the LHC and a future ILC, showing that these exotic Higgs decays can yield sensitivity to unexplored regions of parameter space for ALPs and dark matter models.

Introduction. The Higgs boson discovered at the Large Hadron Collider (LHC) offers a unique window into new physics, and it is paramount to study its properties with precision. Exotic Higgs decays, i.e. decays of the Higgs boson not present in the Standard Model (SM), constitute a primary avenue to probe the existence of new physics [1]. In the last years, there has been an intense experimental program at the LHC to search for such exotic Higgs decays [2-13] (see also [14] and references therein). Such searches have mainly targeted either fully visible final states, e.g. h → 2f 2f (with f, f SM fermions) or a fully invisible Higgs decay (so-called invisible Higgs width). Considering all/part of the Higgs boson decay products in these exotic decays to be invisible at colliders is well-motivated theoretically, e.g. if the Higgs boson directly interacts with a dark (i.e. not feeling the SM gauge interactions) sector of Nature, possibly containing the dark matter (DM) particle(s), or if the Higgs decay products are very long-lived and decay outside the LHC detectors. Yet, partially invisible (semi-dark ) Higgs boson decays constitute a much less explored avenue to search for new physics beyond the SM (BSM) coupled to the Higgs boson, and studies of these semi-dark Higgs decays exist in the literature for very few BSM scenarios [15,16] (see [17][18][19][20][21] for existing experimental searches). Such searches are fully complementary to searches for invisible Higgs decays, and generally probe different regions of parameter space of the same BSM theories. Semi-dark Higgs decays allow in particular to obtain key information on the nature of the coupling between the Higgs and the invisible state(s), by reconstructing the visible part of the exotic Higgs decay.
In this work we target the previously unexplored semidark Higgs decay h → ZX, with X a BSM particle invisible at the LHC (manifesting as missing transverse energy / E T ), and we show it is a promising avenue to probe various well-motivated BSM scenarios: X could be an axion like particle (ALP) or dark photon that decays invisibly or is long-lived, escaping the detector. It could also be a pseudoscalar mediator particle between the SM and a dark sector of Nature containing the DM particle 1 .
We focus our study on the leptonic decay of the Z boson, Z → (with = e, µ), leading to the Higgs final state h → ZX → + / E T . Incidentally, the SM decays h → ZZ * → νν and h → W W * → ν ν yield the same final state. For the latter, the two leptons do not reconstruct the Z boson mass m Z 91 GeV in general, which can be used to tell apart h → ZX from this SM decay process. However, h → ZZ * → νν completely mimics a possible BSM signal. The SM decay h → Zνν then constitutes a "neutrino floor " 2 to experimental searches for new physics in the semi-dark h → ZX (X → / E T ) channel, below which a possible BSM signal would be buried. It also provides a target sensitivity for the h → ZX (X → / E T ) search at the LHC and future colliders which would guarantee a detection (albeit in that case not of BSM physics!), given by the SM branching fraction BR(h → Zνν) SM = 5.4·10 −3 [30].
LHC searches for h → ZX → + / E T . We consider for the rest of this work an LHC center-of-mass energy √ s = 14 TeV. Our analysis reveals the convenience of focusing on the production of the Higgs boson at the LHC in association with a Z boson, pp → Zh: For gluon-fusion (ggF) and vector boson fusion (VBF) Higgs production channels, the Higgs is either produced on its own (ggF) or recoiling against jets (ggF, VBF). Since the phase space for the Higgs decay h → ZX is fairly small (as (m h − m Z )/m h 1), an accurate / E T reconstruction may be limited by the transverse momentum (p T ) 1 A pseudoscalar mediator would nicely explain the absence of a spin-independent signal in current DM direct detection experiments [22]. These pseudoscalar portal to DM scenarios have also been proposed to explain [23][24][25] the γ-ray excess [26,27] in the Fermi-LAT observations of the Milky Way Galactic Center [28]. 2 In analogy to DM direct detection experiments, where coherent elastic neutrino-nucleus scattering can pose an irreducible background to DM searches, known as the "neutrino floor" [29]. resolution of the jets. In addition, the + / E T + jets final state has very large SM backgrounds, in particular reducible ones if the / E T reconstruction is not perfect. Higgs production in association with an electroweak gauge boson, pp → W ± h and pp → Zh, is thus better suited for the h → ZX (X → / E T ) exotic Higgs decay search at the LHC. Yet, the leptonic decay of the W boson in pp → W h adds / E T to the final state, making it challenging to disentangle this contribution from the Higgs boson decay products. In addition, the LHC cross section for the dominant SM background in this case, pp → W ± Z, is very large, O(50) pb. In contrast, for pp → Zh (h → Z + / E T ) the leptonic decay of both Z bosons offers a sharp reconstruction of the two di-lepton resonances together with an accurate / E T measurement, combined with SM backgrounds that can be efficiently suppressed or are very small to begin with, as we discuss in detail below.
For our analysis, we generate the BSM signal specifically using a FeynRules [31] implementation of the two-Higgs-doublet-model plus pseudoscalar singlet (2HDM + a) extension of the SM (see e.g. [32][33][34]), through the decay h → Za (with a invisible). Our results apply to any two-body Higgs decay h → ZX, X → / E T . Both Z bosons from the signal are considered to decay leptonically, Z → . The relevant SM backgrounds are pp → ZZ → 4 (with / E T appearing via mis-measurements and detector effects), pp → ZZZ, W W Z → 4 + 2ν, pp → ttZ, tW Z → 4 + 2ν + jets, and pp → Zh (h → W W * → 2 + 2ν). We generate signal and SM background event samples with MG5 aMC@NLO [35] (using the NNPDF31 nnlo [36] parton distribution functions) with subsequent parton showering and hadroniza-tion via Pythia 8 [37] and detector simulation via Delphes [38], using the detector card designed for High Luminosity (HL)-LHC studies. We normalize the respective cross sections to their next-to-leading-order (NLO) in QCD values, obtained from the literature [39,40] (for the pp → Zh and pp → ZZ processes the normalization is however performed to the NNLO cross section [30,41]; to avoid known issues at NLO in QCD related to real bquark emission [42,43], tW Z is kept at LO with a negligible impact on our analysis). Selected events are required to contain exactly four reconstructed leptons after detector simulation, comprising two pairs of opposite-sign, same-flavor leptons. Events must pass the single, two or three-lepton trigger requirements from the ATLAS 2018 Trigger menu [44]. When multiple di-lepton combinations satisfying the selection requirements exist, the one (with m i the di-lepton invariant masses) is chosen. Extra hadronic activity is vetoed by rejecting events with either b-tagged jets or hard jets with p T > 50 GeV.
Since the Higgs decay is partially invisible, its invariant mass cannot be fully reconstructed, nor can the dilepton pair from its decay be straightforwardly identified. The latter is key to better exploit the kinematic properties of the BSM signal in the analysis. We may identify the di-lepton pair corresponding to the Z boson from the Higgs decay using the transverse mass M T , given , with E T and p T the missing transverse 3-momentum and Z boson transverse 3-momentum, respectively; a complementary approach would be to select the Z boson closest to E T in the azimuthal plane as the one from the Higgs decay. Figure 1 shows the M T (top) and ∆φ(Z, E T ) (azimuthal angle between E T and the 3-momentum of the di-lepton pair, bottom) distributions for the leptonicallydecaying Z boson from the Higgs decay (labelled as Z 1 ) and the Z boson produced in association with the Higgs (labelled as Z 2 ). To optimally exploit the event kinematic information in identifying Z 1 and Z 2 for the BSM signal, we build a neural network (NN) (two hidden layers, 32 nodes each, using rectified linear unit activation for the hidden layers and a sigmoid function for the output) which takes as input M T and ∆φ(Z, E T ) for both di-lepton pairs. The correct and wrong Z i assignments for the NN training are labelled using generator-level information. The NN is trained with a Monte Carlo sample of 20000 signal events (not used in our subsequent analysis) with m X = 1 GeV, using the Adam algorithm for the optimisation. The efficiency obtained for a correct Z 1,2 choice for the signal is 73%, and the NN is then applied in our sensitivity analysis to both the BSM signal (for m X ∈ [1, 32.5] GeV) and the SM backgrounds. The signal cross section (for BR(h → ZX) = 1, BR(X → / E T ) = 1 and m X = 1 GeV) after the ini-  Our h → ZX LHC signal region must target relatively high-p T Zh associated production, with reconstructed Z-boson resonances for the two di-lepton pairs and the Higgs transverse mass M T from the Z 1 di-lepton candidate. Requiring a moderately large amount of E T , demanding Z 1 to be well-aligned with E T in the azimuthal plane and rejecting events with a large rapidity gap between di-lepton pairs also improves the sensitivity of the analysis. The rich event kinematics (four visible objects in the final state plus the missing transverse energy) indicates that a multivariate approach which accesses the full kinematic information of the events could enhance our BSM signal sensitivity. We use another NN (two hidden layers of 256 nodes, same activation functions and optimisation as before) for the discrimination between signal and SM background, with input variables: E T , m 4 (four-lepton invariant mass), m 1 and m 2 , ∆φ(Z 1 , E T ) and ∆φ(Z 2 , E T ), M T (built from Z 1 ), p 1 T and p 2 T (di-lepton transverse momenta), p 1 T , p 2 T , p 3 T , p 4 T (transverse momenta of the four leptons, ordered from higher to lower) and (p 2 T + / E T )/p 1 T . The NN is trained with an unbalanced Monte Carlo set dominated by ZZ → 4 events, precisely to optimise the rejection of this SM background (as it has by far the largest LHC cross section among SM backgrounds). The NN score θ N N for the m X = 1 GeV signal and relevant SM backgrounds is shown in Figure 2 (for X of spin-1, angular correlations in the Z 1 di-lepton pair mildly differ from the spin-0 X case analyzed here, yet our signal sensitivity results would be nearly unchanged).
Our signal region is defined for HL-LHC as θ N N ≥ 0.997. The resulting signal and SM background efficiencies (evaluated on the NN test sets) are 0.12, 1.5 · 10 −4 , 2.8 · 10 −5 , 0.0013, 0.012, 0.0016 and < 9.4 · 10 −4 , respectively for the signal (with m X = 1 GeV), ZZ → 4 , W W Z, ZZ → 2 + 2τ , ZZZ, tW Z and ttZ. We employ the "Asimov estimate" [45] (since O(s/b) is not small) to derive the 2σ exclusion sensitivity on BR(h → ZX) × BR(X → E T ) with the 3 ab −1 of integrated luminosity from HL-LHC, in the range m X ∈ [1, 32.5] GeV (our results may also be directly extrapolated to the m X → 0 limit). We find our NN results to improve by ∼ 30% the sensitivity of a cut-and-count analysis, and come close to probing the Higgs neutrino floor (for We repeat our analysis for an integrated luminosity of 300 fb −1 , with a less stringent signal region cut θ N N ≥ 0.985 to increase the fraction of signal events surviving the selection. The sensitivity results for 300 fb −1 and 3 ab −1 (HL-LHC) are shown in Figure 3.  [46,47] operating at √ s = 250 GeV would be able to probe BR(h → ZX) down to the Higgs neutrino floor by exploiting several advantages over the LHC search discussed in the previous section: (i) Higgstrahlung e + e − → Zh is now the dominant Higgs production mode. (ii) The e + e − collisions at ILC offer a much cleaner environment (largely void of hadronic activity) and the 3-momenta of the incoming particles is known up to radiative and smearing effects, allowing for full missing momentum reconstruction. (iii) The Higgs recoil mass, constructed from the Z-boson recoiling against the Higgs boson (Z 2 ) as

ILC searches for
√ s, provides a straightforward way to correctly identify Z 1,2 for the BSM signal (M reco built out of the Z-boson from the Higgs decay, Z 1 , will not present any resonant structure).
For our analysis, we specifically consider √ s = 250 GeV with 2 ab −1 of integrated luminosity (90% of it evenly split between the two opposite beam helicities) and beam polarizations being 80% for the electrons and 30% for the positrons respectively [47]. Again, we consider the SM Higgs produced in association with a Zboson, e + e − → Zh for our BSM signal. The relevant SM backgrounds are now e + e − → ZZ (→ 4 , 2 + 2τ ), e + e − → W W Z and e + e − → ZZνν (including VBF initiated contributions; Higgs mediated contributions correspond to the SM Higgs neutrino floor, and are not included) Our signal and background event generation is performed as in the previous section, using in this case the Delphes detector card designed for ILC studies [38,48] (a study of lepton collider capabilities including the effects of initial state radiation or beamstrahlung is left for future work).
Our initial event selection mimics that of the previous LHC analysis. The cross sections for the signal (for BR(h → ZX) = 1 and m X = 1 GeV) and SM backgrounds after event selection are respectively 1.421 fb, 5.64 fb (for ZZ → 4 ), 0.13 fb (for ZZ → 2 2τ ), 0.073 fb (for W W ( * ) Z → 4 + 2ν, dominated by the e − e + → Zh, h → W W * contribution), and 0.011 ab (for ZZνν → 4 + 2ν). For the ILC environment, the use of a NN does not offer such a strong advantage over a simpler (cut-and-count) analysis. We thus define our kinematic region for signal extraction in the latter way: we demand reconstructed Z-boson resonances for the two dilepton pairs, m Z1 ∈ [55, 100] GeV, m Z2 ∈ [80, 105] GeV; we require the recoil mass constructed out of Z 2 in the range M reco ∈ [120, 135] GeV, together with p Z2 > 50 GeV and / E ∈ [5, 60] GeV (p Z2 , / E respectively the modulus of the Z 2 di-lepton candidate's 3-momentum and the modulus of the missing 3-momentum); the invariant mass m miss Z1 built from Z 1 and / E, given by (m miss GeV) 2 . These signal region cuts have an efficiency of 0.89 for the BSM signal (for m X = 1 GeV), and 1.7 · 10 −5 , 0.013, 0.085, 0.24 for the ZZ → 4 , ZZ → 2 2τ , W W Z and ZZνν SM backgrounds, respectively. Using the "Asimov estimate", we derive a 2σ sensitivity BR(h → ZX) × BR(X → / E) = 0.0045. We show the corresponding ILC sensitivity as a function of m X in Figure 3. We note that this sensitivity lies below the SM Higgs neutrino floor, not been included in our analysis. This means that we should now instead consider the ILC discovery potential of BR(h → Zνν) SM ; the expected significance over the background only hypothesis reaches 2.4σ (the significance may be enhanced to ∼ 4σ, at the expense of our BSM analysis not yielding a uniform sensitivity in m X ). ILC can thus sweep the entire new physics parameter space of semi-dark Higgs decays down to the Higgs neutrino floor.
Constraints on specific models. (i) Axion-like particles. The existence of interactions between the SM Higgs and light axion-like particles (ALPs) is a well-motivated BSM possibility [49,50]. Exotic Higgs decays represent a key experimental signature in this case. The Higgs boson partial decay width into a SM Z-boson and the ALP a is Γ(h → Za) = (m 3 h /16πf 2 a ) c 2 aZh λ 3/2 , with f a the ALP decay constant, a /m 4 h and c aZh the Wilson coefficient for the effective operator that couples the ALP to the SM Higgs field (which we leave here unspecified, see [49,50] for details). If a then couples to some hidden sector particle(s) (see e.g. [51,52]), its dominant decay mode(s) may be into the dark sector, thus invisible at colliders. This encompasses the intriguing possibility that the ALP may be a mediator between the SM and the DM candidate [51]. Higgs decays h → Za, a → E T can then probe such ALP scenarios. To translate the model-independent LHC and ILC projected sensitivities from the previous sections into a probe of the parameter space of ALPs, we specifically consider, together with the coupling between the ALP and the SM Higgs, an ALP coupling to a hidden fermion χ, given by y χχ γ µ γ 5 χ ∂ µ a/f a as well as an ALP coupling to photons c aγγ /f a a F µνF µν (we do not include an ALP coupling to gluons or SM fermions for simplicity). y χ does not have a preferred value, while the expectation for the bosonic Wilson coefficient is c aγγ ∼ α EM (the electromagnetic coupling constant) [53]. We then set c aγγ = α EM (Q) (Q being the energy scale of the process considered), and y χ = 1, c aZh = 1, m χ = 0.45 m a (to allow for the invisi-ble ALP decay a → χχ), and show in Figure 4 the LHC projected 2σ sensitivity on BR(h → Za) × BR(a → χχ) in the (m a , f a ) plane. We also depict the Higgs neutrino floor, within reach of the √ s = 250 GeV ILC. Figure 4  also shows, under the assumption that χ is the DM particle, the (m a , f a ) relation yielding (for the choice of parameters described above) the observed DM relic abundance Ω DM h 2 = 0.12 [54] generated via thermal freezeout in the early Universe (taken from [51]), as well as the existing and projected constraints on this ALP scenario from searches at LEP, LHC and flavor factories (Babar, Belle-II), and astrophysics (supernova 1987A), all detailed in Appendix I. Finally, we also show in Appendix I the corresponding limits on the (m a , f a ) plane if a hypercharge coupling c aBB /f a a B µνB µν (rather than a coupling only to photons) is assumed for the ALP.
(ii) 2HDM + a. Two Higgs doublet models (2HDM) extended by a singlet pseudoscalar mediator and a fermionic singlet DM particle constitute the minimal renormalizable realization of a pseudoscalar portal to DM [25,[32][33][34]55], which avoids the stringent existing DM direct detection constraints [22] (it yields a spinindependent DM-nucleon scattering cross section only at loop level [56,57]); it can also fit the observed gammaray galactic center excess [25,58]; and it is a leading benchmark scenario for the DM interpretation of LHC searches [59]). The potential of the 2HDM+ a is [34] with real pseudoscalar mediator a 0 and Dirac fermion DM χ with mass m χ , both singlets under the SM gauge interactions. V 2HDM is the 2HDM scalar potential (with a Z 2 -symmetry softly-broken by a µ 2 12 H † 1 H 2 +h.c. term in V 2HDM , see e.g. [60]). The κ term in (1) yields a mixing θ between a 0 and the 2HDM neutral pseudoscalar state A 0 , resulting in two mass eigenstates a, A (with m a < m A ) which provide the portal between the SM and DM.
For a light pseudoscalar a (singlet-like), the coupling between a, h and Z leads to semi-dark Higgs decays (a decays to DM particles with a branching fraction BR(a → χχ) 1 unless y χ 1). The partial decay width is Γ(h → Za) = (m 3 h /16πv 2 ) c 2 β−α sin 2 θ λ 3/2 , with v the electroweak vev and c β−α ≡ cos(β − α) parametrizing the deviation from the 2HDM alignment limit [60] (with tan β = v 2 /v 1 , vev ratio of 2HDM Higgs doublets, and α the mixing angle between the 2HDM CP-even weak eigenstates [61,62]). The model also features an h → aa decay, leading to a Higgs invisible partial width via a → χχ decays. For |c β−α | 1, as needed to satisfy the present LHC Higgs signal strength measurements [63,64], we generally expect Γ(h → aa) Γ(h → Za), yet in certain (albeit small) regions of the 2HDM+ a parameter space, the h → Za semi-dark Higgs decay can provide stronger sensitivity than the h → aa invisible Higgs decay, in particular when the h − a − a coupling (see Appendix II for details) vanishes. We consider here a Type-I 2HDM (see [62]) with c β−α = 0.2, t β ≡ tan β = 6, M = 600 GeV, m H0 = m H ± = m A0 = 700 GeV (with M 2 = µ 2 12 /s β c β and with H ± and H 0 the charged and neutral CP-even heavy 2HDM scalars). We choose in addition m χ = 0.45 m a and y χ fixed to yield the observed DM relic density via thermal freeze-out (see e.g. [33,58] and Appendix II). We further consider λ a1 = λ a2 fixed in each case to the value that yields Γ(h → aa) = 0, and show in Figure 5 the projected LHC sensitivity (with 300 fb −1 and at HL-LHC with 3 ab −1 ) of the semi-dark Higgs decay h → Z( )+ E T in the (m a , sin θ) parameter space of the 2HDM+ a. In addition, while our 2HDM (c β−α , t β ) benchmark satisfies both present and HL-LHC projected limits from Higgs signal strengths on 2HDM parameters [65], we find using the ScannerS [66] and Hig-gsSignals [67,68] numerical codes that h → Za exotic Higgs decays in this scenario are currently constrained at 95% C.L. to BR(h → Za) < 0.042 from Higgs signal strength measurements. This limit is also shown in Figure 5. Further details on these constraints, as well as a discussion on other searches that could be sensitive to the 2HDM+ a parameter space (yet not to the specific benchmark we choose here) are given in Appendix II. (iii) A comment on dark photons. Light dark photons Z D [69] which interact with the SM via kinetic mixing (see e.g. [70]) give rise to an exotic Higgs decay h → ZZ D . Invisible dark photons would then constitute another new physics scenario that semi-dark Higgs decays could be sensitive to. However, current 95% C.L. bounds on the kinetic mixing parameter from EW precision observables set < 0.03 for dark photon masses < 30 GeV [71]. The corresponding h → ZZ D branching fraction is then < 10 −3 (see e.g. [1]), below the Higgs neutrino floor.  ITN HIDDeN), as well as from the AEI through the grant IFT Centro de Excelencia Severo Ochoa CEX2020-001007-S.

Appendix I: Details on ALP constraints
We here discuss the constraints on the (m a , f a ) parameter space of an invisibly-decaying ALP under the assumptions made in the main text. From [51], we obtain the 95% C.L. LEP limits from mono-photon searches [72,73], which constrain an ALP produced via its coupling to photons and decaying invisibly, as well as the 90% C.L. limits from e + e − → γ + / E (with / E the missing energy of the event) and rare upsilon decays into γ + / E from Babar [74,75] (see also [76]) and the projected 90% C.L. sensitivity of Belle-II in the γ + / E final state. Also shown in Figure 4 are the current 95% C.L. limits from heavy-ion (Pb-Pb) collisions at the LHC, from ALP searches in light-by-light scattering (as proposed in [77]) performed by ATLAS [78] and CMS [79] (see also [80]). We also include a projection drawn from rescaling the current ATLAS expected sensitivity to an (optimistic) integrated luminosity of 30 nb −1 . Finally, we show the bound from the energy loss of supernova 1987A from ALP emission, as taken from [51]. The supernova 1987A limit is stronger than usually quoted for an ALP coupling only to photons since the invisible decay of the ALP allows its corresponding energy to escape the supernova core even for parameter regions with a sizable coupling to photons (contrary to the usual case, where a large enough coupling to photons will result in the ALP being trapped in the core [81]).
We note that the existence of the invisible decay mode of the ALP leads (under the assumptions discussed in the main text) to a strongly suppressed ALP branching fraction to two photons, BR(a → γγ) ∼ 3 × 10 −4 . Limits from ALP searches in visible final states, like triphoton searches at LEP 1 and LEP 2 via the process e + e − → γ * → γ a, a → γγ (as studied in [82,83]), and searches for Z → γγ decays 3 at LEP 1 have to be rescaled by BR(a → γγ) (assumed to be 100% in [82,83]), and are too weak to appear in Figure 4. Similarly, the dominant invisible decay of the ALP significantly weakens the limits from beam-dump experiments as compared to the case of visible ALP scenarios (see e.g. [84]) roughly by a factor BR(a → γγ) ∼ 10 −4 . A naive re-scaling of beam-dump limits results in no meaningful constraints (beyond what is currently excluded by other experiments/observations) from these, and we choose not to include them in Figure 4. A precise re-derivation of these limits requires to additionally take into account the geometry of each experiment to obtain the expected number of a → γγ events in the detector decay volume, which is beyond the scope of this paper. For an ALP coupled to the hypercharge field strength via c aBB /f a a B µνB µν , this introduces a coupling of the ALP to ZZ and Zγ (besides the already considered coupling to photons), given respectively by c aZZ /f a a Z µνZ µν and c aZγ /f a a F µνZ µν , with c aZZ = sin 2 θ W c aBB and c aZγ = −2 sin θ W cos θ W c aBB (with θ W the weak mixing angle). Fixing c aBB = α EM /cos 2 θ W to match the ALP coupling to photons c aγγ we assume in the main text, the above limits do not change, yet from c aZγ = 0 there is another constraint from LEP searches for rare Z → γ + a decays, with a invisible. The L3 collaboration at LEP has set a limit BR(Z → γ + a) < 1.1 × 10 −6 at 90% C.L. [85], shown in Figure 6 (in purple) together with the already considered constraints on our scenario (Figure 4). Other potential bounds from rare Z decays at LEP and LHC, e.g. from Z → 3γ or Z → γ (see [50]), do not lead to meaningful constraints in Figure 6.
Finally, we comment on the possibility of probing the ALP a via exotic Higgs decays h → Za with a → γγ, as discussed in [50]. We note that, while the corresponding final state allows to consider Higgs production in gluonfusion (resulting in an O(50) enhancement of the Higgs production cross section w.r.t. our scenario, which must rely on Higgs associated production), this is counteracted by the large suppression from BR(a → γγ). A preliminary study of the LHC sensitivity to ALPs in exotic Higgs decays via pp → h → Za, Z → , a → γγ including the leading SM backgrounds has been performed in [86], indicating that such decay channel is much less sensitive than the one discussed in this work (given our assumptions for the ALP branching fractions).

Appendix II: Details on 2HDM + a constraints
The main direct experimental probes of the existence of a light (m a < 30 GeV) singlet-like pseudoscalar a in the 2HDM + a are the exotic Higgs decays h → Za (which this work explores in detail) and h → aa. For the latter, the partial width is Γ(h → aa) = (v 2 /32 π m h ) × g 2 haa 1 − 4m 2 a /m 2 h , with the h − a − a coupling g haa given in Eq. (2). In the main text we have analyzed the tuned scenario g haa = 0 (with varying λ a1 = λ a2 ). We now consider the same 2HDM parameter benchmark as in the main text (c β−α = 0.2, t β = 6, M = 600 GeV, m H0 = m H ± = m A0 = 700 GeV), but fix λ a1 = λ a2 = 0 such that Γ(h → aa) = 0. The resulting modified LHC sensitivity to the (m a , sin θ) parameter space of the 2HDM+ a from semi-dark Higgs decay h → Z( ) + E T is shown in Figure 7. We also depict in Figure 7 the constraint on the Higgs invisible width from h → aa decays, which at present is BR(h → E T ) < 0.11 at 95% C.L. [87] under the assumption of SM Higgs production, and is expected to be BR(h → E T ) < 0.04 at 95% C.L. [88] at the HL-LHC.
Other searches for the state a do not provide meaningful sensitivity in the scenario we consider: searches for h → aa in visible final states (see [14] for a review) like bbτ τ [10] and τ τ τ τ [6] are found to be O(10 3 ) less sensitive than probes of the Higgs invisible width, and fall short of providing any limit on BR(h → aa) by a factor ∼ 50 − 100, with searches in other final states (e.g. bbµµ [8], τ τ µµ [4, 5]) yielding even smaller sensitivity. Such visible decays of a are then generally not relevant in the 2HDM + a with m a > 2 m χ , since matching the observed DM relic density requires y χ ∼ 1 (see below), leading to BR(a → χχ) > 0.99 in general. We also find that current LHC mono-jet searches [89] fall short of probing any region of the (m a , sin θ) plane of Figures 5 and 7 by a factor ∼ 1/(3 t 2 β ). Finally, we note that, while the a − h − Z coupling could be constrained via Higgs boson production in association with missing energy at LEP2 through the process e + e − → Z * → ah, this does not yield meaningful limits on the 2HDM + a parameter space: the searches for h +νν signatures by the OPAL [90], L3 [91] and ALEPH [92,93] experiments at LEP impose a constraint on the missing mass M miss of the event (equal to m a in our scenario) which is not fulfilled by our signal (e.g. the OPAL analysis [90] requires 50 GeV < M miss < 130 GeV, and it is thus not sensitive to m a < 30 GeV). The corresponding search by the DELPHI experiment [94], while not imposing such a cut on M miss , does not consider Higgs boson masses above 120 GeV.
We also consider possible constraints on the spin-0 states from the 2HDM, H ± , H 0 and A (doublet-like). Electroweak precision observables (EWPO) constrain (dominantly via the oblique T -parameter) the mass splittings among the 2HDM scalars, since the 2HDM scalar potential of the 2HDM is custodially invariant for m A0 = m H ± or m H0 = m H ± . The latter is chosen for our benchmark scenario, directly satisfying EWPO. Potential constraints from flavour physics, in particular from B-meson decays:B → X s γ decays (which constrain m H ± [95,96]) and from B s → µ + µ − (which constrain the presence of a light a coupling to SM fermions [97,98]), are also directly satisfied for a Type-I 2HDM with moderately high t β (we have chosen in particular t β = 6 in our benchmark analysis). Finally, we also discuss direct searches for the 2HDM states as a probe of our scenario: away from the 2HDM alignment limit, H 0 → W + W − and H 0 → ZZ decays could yield sensitivity if the mass scale of the 2HDM scalars is not very high. At the same time, H 0 → Za decays would lead to resonant mono-Z signatures [32,34], which have been recently been searched for by the AT-LAS experiment [99] (ATLAS also searches for this final state in H 0 → ZZ → νν decays [100]). However, in all these cases, we find that for |s β−α t −1 β − c β−α | 1 (as our scenario features) the production of H 0 at the LHC is suppressed, and no meaningful limit is obtained. To conclude, we discuss the need to reproduce the observed DM relic density Ω DM h 2 = 0.12 [54] within the 2HDM + a, via DM thermal freeze-out in the early-Universe. For m χ 2 GeV, the DM annihilation cross section into SM particles (via χχ → q SM q SM , with q SM here being generic SM particles) in the nonrelativistic limit is with Γ a the decay width of a. The sum over SM fermion f annihilation channels involves quarks (N C = 3) and charged leptons (N C = 1). Reproducing the observed DM relic density via thermal freeze-out requires σv 3 × 10 −26 cm 3 /s, which generically leads to O(1) values for y χ . For m χ < 2 GeV the DM annihilation into SM fermions (b, c-quarks and/or τ -leptons, depending on m χ ) ceases to be the dominant DM annihilation process, and instead annihilation into QCD hadrons (via the 1-loop coupling of a to gluons) dominates. Due to its complexity, we do not study that region of the 2HDM + a parameter space, which may also involve y χ 1, in this work. * ja.