Mixed hidden sector/visible sector dark matter and observation of CP odd Higgs at HL-LHC and HE-LHC

It is very likely that similar to the case of visible matter, dark matter too is composed of more than one stable component. In this work we investigate a two-component dark matter with one component from the visible sector and the other from the hidden sector. Specifically we consider a $U(1)_X$ hidden sector extension of MSSM/SUGRA where we allow for kinetic and Stueckelberg mass mixing between the two abelian $U(1)'s$, i.e., $U(1)_X$ and $U(1)_Y$. We further assume that the hidden sector has chiral matter which leads to a Dirac fermion as a candidate for dark matter. The lightest neutralino in the visible sector and the Dirac fermion in the hidden sector then constitute the two components of dark matter. We investigate in particular MSSM/SUGRA models with radiative breaking occurring on the hyperbolic branch where the Higgs mixing parameter $\mu$ is small (order the electroweak scale) which leads to a lightest neutralino being dominantly a higgsino. While dark matter constituted only of higgsinos is significantly constrained by data on dark matter relic density and by limits on spin independent proton-DM scattering cross section, consistency with data can be achieved if only a fraction of the dark matter relic density is constituted of higgsinos with the rest coming from the hidden sector. An aspect of the proposed model is the prediction of a relatively light CP odd Higgs $A$ (as well as a CP even $H$ and a charged Higgs $H^{\pm}$) which is observable at HL-LHC and HE-LHC. We perform a detailed collider analysis search for the CP odd Higgs using boosted decision trees in $\tau_h\tau_h$ final states and compare the discovery potential at HL-LHC and HE-LHC. We show that while several of the points among our benchmarks may be observable at HL-LHC, all of them are visible at HE-LHC with much lower integrated luminosities thus reducing significantly the run time for discovery.


Introduction
The discovery of the Higgs boson at the Large Hadron Collider (LHC) at ∼ 125 GeV [1,2], gives strong support for supersymmetry (SUSY). This is so because within the Standard Model (SM) the Higgs boson mass can lie in a wide range up to several hundred GeV in mass [3], while within supersymmetry/supergravity (SUGRA) unified models (for a review see, e.g. [4]) the mass of the Higgs boson is predicted to lie below 130 GeV [5][6][7]. In addition to the fact that LHC data respects the supersymmetric limit, stability of the vacuum can be preserved within supersymmetry up the Planck scale while within the standard model the vacuum stability holds only till around 10 10 GeV [8,9]. However, as is well known the Higgs boson mass at ∼ 125 GeV requires a large loop correction within the minimal supersymmetric standard model MSSM/SUGRA which in turn implies that the size of weak scale supersymmetry is large, lying in the several TeV region. This also explains why SUSY has not been observed at accelerators thus far. The large size of weak scale supersymmetry also has implications for dark matter (DM). In SUGRA unified models, the low energy sparticle spectrum determined by running the renormalization group equations (RGE) and radiative breaking of the electroweak symmetry has several branches. On one branch, the ellipsoidal branch, the Higgs mixing parameter µ is large and the lightest supersymmetric particle (LSP), the lightest neutralino, is typically a bino. However, in the early universe the binos are not annihilated efficiently leading to an LSP relic density far in excess of the current experiment [10]. Here one way to reduce the relic density is through the utilization of coannihilation. Other possibilities within SUGRA models to get conformity with the relic density constraint include a wino-like dark matter or a higgsino-like dark matter.
In this work we focus on SUGRA models on the hyperbolic branch with a small µ (of order the electroweak scale) where the LSP is higgsino-like. Indeed within radiative breaking of the electroweak symmetry higgsino-like dark matter can arise naturally on the hyperbolic branch when µ is small [11][12][13] (for related works, see, e.g., [14][15][16][17]). Models of this type are severely constrained by simultaneous satisfaction of dark matter relic density data and by the spin-independent proton-DM scattering cross section limits in direct detection experiments. However, such models can be viable if dark matter is multicomponent with the higgsino-like DM contributing only a fraction of the relic density with the remainder made up from other sources. Here we discuss a two-component dark matter model where one component is the higgsino (a Majorana fermion) of the visible sector while the other component arises from the hidden sector and is a Dirac fermion [18]. Thus the two-component dark matter model is a U (1) X extension of the standard model gauge group where the U (1) X gauge boson of the hidden sector and the U (1) Y gauge boson of the visible sector have both kinetic [19,20] and Stueckelberg mass mixings [21][22][23] (for Stueckelberg extension with an enlarged gauge group, see [24,25]). Further, the hidden sector contains matter which provides a Dirac fermion as the second component of dark matter. It is then seen that the Dirac fermion of the hidden sector provides the dominant piece of the relic density but the higgsino dark matter dominates the spin independent cross section in the direct detection experiments.
One remarkable aspect of the two-component model is the prediction of a relatively light CP odd Higgs (in the range of few hundred GeV) which lies in the observable range of the future generation of colliders (see, e.g., [26,27] and [28,29]). Specifically we focus here on the high luminosity LHC (HL-LHC) and high energy LHC (HE-LHC). In this work we carry out a detailed analysis of the integrated luminosities needed for the observation of this low-lying Higgs. Its observation would lend support to the higgsino nature of the LSP and the multi-component nature of dark matter. At the same time some of the predicted spinindependent scattering cross-sections also lie in the range of the next generation dark matter direct detection experiments. The outline of the rest of the paper is as follows: Details of the two-component model are discussed in section 2. The scalar sector of the theory is further elaborated in section 3. In section 4 we give ten representative benchmarks satisfying the relic density constraint along with the Higgs boson mass constraint. An analysis of the two-component dark matter, of relic density and of direct detection is discussed in section 5. Associated production of CP odd Higgs along with heavy quarks at the LHC is discussed in section 6 followed by the prospects of discovering a CP odd Higgs at HL-LHC and HE-LHC in section 7. Conclusions are given in section 8. It is also shown that part of the parameter space of the extended model can be probed in the next generation direct detection experiments such as XENONnT and LUX-ZEPLIN.
We note in passing that there are a variety of supersymmetric U (1) extensions and their effect on DM and collider analyses have been studied extensively in the literature [30,31]. We also note that the two-component model can be easily extended to include other forms of dark matter such as an axion [32] or an ultralight axion [33][34][35]. Further, several works on the HL-LHC and HE-LHC discovery potential have appeared recently and in the past years [36][37][38][39][40][41]].

The model
As discussed above we consider an extension of the standard model gauge group by an additional abelian gauge group U (1) X of gauge coupling strength g X . The MSSM particle spectrum in the visible sector, i.e., quarks, leptons, Higgs and their superpartners are assumed neutral under U (1) X . Thus the abelian gauge sector of the extended model contains two vector superfields, a vector superfield B associated with the hypercharge gauge group U (1) Y , a vector superfield C associated with the hidden sector gauge group U (1) X , and a chiral scalar superfield S. In the Wess-Zumino gauge the B and C superfields have the following components and The chiral scalar superfield S can be expanded in terms of its component fields as The gauge kinetic energy sector of the model is given by Next we allow gauge kinetic mixing between the U (1) X and U (1) Y sectors with terms of the form As a result of Eq. (5) the hidden sector interacts with the MSSM fields via the small kinetic mixing parameter δ. The kinetic terms in Eq. (4) and Eq. (5) can be diagonalized by the transformation Aside from gauge kinetic mixing, we assume a Stueckelberg mass mixing between the U (1) X and U (1) Y sectors so that [21] We note that Eq. (7) is invariant under U (1) Y and U (1) X gauge transformation so that, In component notation, L St is In the unitary gauge, the axion field a is absorbed to generate mass for the U (1) X gauge boson.
The matter sector of the model consists of the visible sector chiral superfields denoted by Φ i where i runs over all quarks, squarks, leptons, sleptons, Higgs and Higgsino fields of the MSSM and hidden sector chiral superfields denoted by Ψ i . The Lagrangian for the matter interacting with the U (1) gauge fields is given by where Y is the U (1) Y hypercharge and X is the U (1) X charge. The MSSM fields are not charged under the hidden sector and vice-versa, i.e., XΦ i = 0 and Y Ψ i = 0. The minimal particle content of the hidden sector consists of a left chiral multiplet Ψ = (φ, f, F ) and a charge conjugate Ψ c = (φ , f , F ) so that Ψ and Ψ c carry opposite U (1) X charge and hence constitute an anomaly-free pair. The Dirac field ψ formed by f and f has a mass M ψ arising from the term M ψ ΨΨ c in the superpotential. Following SUSY breaking, the scalar fields of the hidden sector acquire soft masses equal to m 0 (the universal scalar mass of the visible sector) and the additional Dirac mass such that It is convenient from this point on to introduce Majorana spinors ψ S , λ X and λ Y so that In addition to the MSSM soft SUSY breaking terms, we add new terms pertinent to the additional fields where m X is the U (1) X gaugino mass and M XY is the U (1) X − U (1) Y mixing mass.
After electroweak symmetry breaking, ψ S and λ X mix with the MSSM gauginos and higgsinos to form a 6 × 6 neutralino mass matrix. We choose as basis (λ Y , λ 3 ,h 1 ,h 2 , λ X , ψ S ) where the last two fields arise from the extended sector and the first four, i.e., λ Y , λ 3 ,h 1 ,h 2 are the gaugino and the higgsino fields of the MSSM sector. Using Eq. (6) we rotate into the new basis (λ Y , λ 3 ,h 1 ,h 2 , λ X , ψ S ) so that the 6 × 6 neutralino mass matrix takes the form where s β ≡ sin β, c β ≡ cos β, s W ≡ sin θ W , c W ≡ cos θ W with M Z being the Z boson mass. We label the mass eigenstates asχ whereχ 0 5 andχ 0 6 belong to the hidden sector and mix with the usual MSSM neutralinos. In the limit of small mixings between the hidden and the MSSM sectors the masses of the hidden sector neutralinos are We turn now to the charge neutral gauge vector boson sector. Here the 2 × 2 mass square matrix of the standard model is enlarged to become a 3×3 mass square matrix in the U (1) Xextended SUGRA model. Thus after spontaneous electroweak symmetry breaking and the Stueckelberg mass growth the 3 × 3 mass squared matrix of neutral vector bosons in the basis (C µ , B µ , A 3 µ ) is given by where A 3 µ is the third isospin component, The mass-squared matrix of Eq. (17) has one zero eigenvalue which is the photon while the other two eigenvalues are where M + is identified as the Z boson mass while M − as the Z boson. The diagonalization of the mass-squared matrix of Eq. (17) can be done via two orthogonal transformations where the first is given by [23] which transforms the mass matrix to M 2 where = c δ − s δ . The gauge eigenstates of M 2 V can be rotated into the corresponding mass eigenstates (Z , Z, γ) using the second transformation via the rotation matrix with c η (c θ )(c φ ) ≡ cos η(cos θ)(cos φ) and s η (s θ )(s φ ) ≡ sin η(sin θ)(sin φ), where η represents the mixing angle between the new gauge sector and the standard model gauge bosons while the other angles are given by . The resulting mixing angle is thus given by 3 The scalar sector of the U (1) X -extended MSSM/SUGRA The addition of the chiral scalar superfield S and the hidden sector matter fields bring about new scalar fields to the theory. Thus, the scalar fields of the U (1) X -extended MSSM/SUGRA are the Higgs fields, the scalar ρ and the fields φ and φ of the hidden sector. In the MSSM, the Higgs sector contains two Higgs doublets H d and H u , with opposite hypercharge which ensures the cancellation of chiral anomalies. Here H d gives mass to the down-type quarks and the leptons while H u gives mass to up-type quarks. The Higgs potential in the MSSM arises from three sources: the F term of the superpotential, the D terms containing the quartic Higgs interaction and the soft SUSY breaking Higgs mass squared, m 2 H d and m 2 Hu , and the bilinear B term. The additional scalar field ρ enters the Higgs potential and mixes with the MSSM Higgs doublets. The full CP-conserving Higgs scalar potential in the extended model can be written as where µ is the Higgs mixing parameter appearing in the superpotential term µĤ u ·Ĥ d . The neutral components of the Higgs doublets and the scalar ρ can be expanded around their VEVs so that The MSSM 2 × 2 Higgs mass matrix is now extended to become 3 × 3 with the new scalar field φ ρ mixing with the two CP even Higgs fields φ d and φ u . As a result, the masses of the CP even Higgses h and H are modified by amounts proportional to and δ which, however, are small. Similarly, the corrections to the CP odd Higgs mass induced by the new sector are negligible. Minimizing the Higgs potential of Eq. (25) in the φ d , φ u and φ ρ directions we obtain the constraints where The last of Eqs. (27) gives which is typically small since and δ are small.

U (1) X -extended MSSM/SUGRA benchmarks
The particle content of the U (1) X -extended MSSM/SUGRA discussed in sections 2 and 3 consists of the particles of the MSSM, and from the hidden sector three spin zero particles (ρ, φ, φ ), three spin 1/2 particles (a Dirac fermion ψ and two Majorana neutralinosχ 0 5 ,χ 0 6 ) and one massive vector boson Z . The model is implemented in the Mathematica package SARAH v4.14.1 [42,43] which generates model files for SPheno-4.0.3 [44,45] which in turn produces the sparticle spectrum and CalcHep/CompHep [46,47] files used by micrOMEGAs-5.0.4 [48] to determine the dark matter relic density and UFO files which are input to MadGraph5 [49]. The input parameters of the U (1) X -extended MSSM/SUGRA [50] with hidden sector matter are taken to be m 0 , A 0 , m 1 , m 2 , m 3 , M 1 , m X , M ψ , B ψ , δ, g X , tan β, sgn(µ), where m 0 , A 0 , m 1 , m 2 , m 3 , tan β and sgn(µ) are the universal scalar mass, the trilinear coupling, the U (1), SU (2) and SU (3) gaugino masses of the MSSM sector and B ψ is the bilinear parameter of the Dirac mass term in the superpotential all taken to be at the GUT scale. Table 1 shows ten representative benchmarks covering a mass range of the CP odd Higgs from ∼ 300 GeV to 750 GeV. It is to be noted that the larger fraction of the relic density is contributed by the Dirac fermion of the hidden sector while the higgsino-like neutralinos contribute smaller fraction. The CP odd Higgs mass along with the neutralino, chargino, stop, gluino and stau masses are presented in Table 2. We also show the mass of the lightest CP even Higgs consistent with the observed 125 GeV Higgs within ±2 GeV error. In some of those benchmarks, the value of m 0 is quite small, for instance, point (g) has m 0 ∼ 800 GeV while the stop mass is ∼ 5 TeV. The reason is the large value of m 3 which via the RGE running generates squark masses in the several TeV range [51]. With heavy gluinos and stops, experimental limits on their masses from ATLAS and CMS can be evaded. Also, the LSP and chargino masses presented in Table 2 have not yet been ruled out by experiment.  Table 2: Display of the SM-like Higgs boson mass, the stau mass, the relevant electroweak gaugino masses, the CP odd Higgs mass and the relic density for the benchmarks of Table 1 computed at the electroweak scale. All masses are in GeV.
Similar MSSM benchmark scenarios have appeared in [52] where heavy Higgses can be probed at the LHC using an effective field theory approach with a spectrum containing light charginos and neutralinos while the rest of the SUSY particles are heavy.
The particle spectrum of the model contains an extra neutral massive gauge boson, Z . Stringent constraints are set on the mass of the Z [53] and most recently by ATLAS [54] using 139 fb −1 of data. In new physics models containing Z with SM couplings, the mass limit is set at m Z > 5.1 TeV. For a model with an extra U (1) X with a gauge coupling strength g X , the limit can be written as For the benchmarks of Table 1, the Z mass obtained from Eq. (18) is ∼ M 1 since M 2 ∼ 0 and s δ 1. Thus the spectrum contains a Z with a mass range of ∼ 800 GeV to ∼ 1500 GeV. However, since the U (1) X coupling g X < 0.1, the limit of Eq. (30) is satisfied for all the benchmarks. The smallness of g X also means that the Z coupling to SM particles is tiny, therefore its production cross-section at pp colliders is suppressed and thus a Z in the mass range noted above is consistent with the current experimental constraints. Further, from Eq. (18), the Z boson mass receives a correction due to gauge kinetic and mass mixings.

Two-component dark matter and its direct detection
As noted earlier one of the constraints on dark matter models is the relic density constraint which according to the PLANCK collaboration [10] is given by Since there are two components to dark matter, the total relic density is the sum of the neutralino and the Dirac fermion relic densities, i.e., Further, the spin independent DM-proton cross-section that enters in direct detection experiments is given by where R χ = (Ωh 2 ) χ /(Ωh 2 ) PLANCK , and R ψ = (Ωh 2 ) ψ /(Ωh 2 ) PLANCK .
Thus one finds that not only the sum of the neutralino and the Dirac fermion relic densities but also their individual contributions have observable consequences as seen from Eq. (34).
The main processes that enter in the neutralino and Dirac fermion relic abundance are Note that the process χψ ↔ SM SM cannot happen since the only allowed vertex is ψχφ and φ does not couple to SM particles. In order to calculate the relic density in the two-component model, one must solve the Boltzmann equations for the χ and ψ number densities n χ and n ψ . Taking into consideration the processes in Eq. (36), the coupled Boltzmann equations for n χ and n ψ are [18] dn where σv χχ denotes σv χχ→SM SM , σv ψψ refers to ψψ → SM SM, χχ and n eq represents the equilibrium number density. The factor of 1/2 appearing in some of the terms is due to the fact that ψ is a Dirac fermion. With the exception of point (e) in Table 2, M ψ > Mχ0 1 and so ψ freezes-out earlier, i.e. at a higher temperature, T f , thanχ 0 1 . Notice the small mass gap betweenχ 0 1 andχ ± 1 resulting in the activation of the coannihilation channel. The solution to Eqs. (37) does not have a closed form and must be solved numerically. The total relic density can be expressed as where x f = m/T f (m being the mass of the DM particle) and the C's are constants proportional to g * −1/2 M −1 Pl with g * the effective number of degrees of freedom at freeze-out for χ or ψ and M Pl is the Planck mass. In columns 2−4 of Table 3 we give the size of the thermally averaged annihilation cross-sections of the processes in Eq. (36). The largest cross-section is that of χχ ↔ SM SM since it has weak scale couplings while reactions involving ψ have cross-sections that are two to five orders of magnitude less than those of χχ. The reason is that ψψ annihilation proceeds through an s-channel exchange of γ, Z, Z with couplings proportional to g X , s δ and . The smallness of those parameters renders the cross-section tiny.  and R ψ × σ SI pψ for the benchmarks of Table 1.
We turn now to the spin-independent (SI) proton-DM scattering cross-section in direct detection experiments. The main contribution to the proton-neutralino scattering cross-section comes from a t-channel exchange of a Higgs boson (h/H) while the proton-Dirac fermion scattering cross-section involves a box diagram with the exchange of the scalar φ from the hidden sector which explains its small value (last column in Table 3) compared to the protonneutralino one (fifth column in Table 3). To get an idea on where the benchmarks lie relative to the current experimental sensitivity of direct detection experiments, we plot in Fig. 1 the ten benchmarks along with the most recent experimental limits from XENON1T [55] for the proton-neutralino cross-section (the color code for those points has no meaning here). We note that some of these benchmarks may be accessible in future experiments such as XENONnT and LUX-ZEPLIN [56]. Note that for most of the benchmarks R ψ × σ SI pψ is very small and lies well below the sensitivity limit of future experiments while others lie even below the neutrino floor.  Figure 1: The SI proton-neutralino cross section exclusion limits as a function of the LSP mass from XENON1T (taken from [55]). The ten benchmarks are overlaid on the plot showing them lying below the upper limit (black curve). The inset shows the limits from LUX 2017 [57], PandaX-II [58] and XENON1T along with the uncertainty bands normalized to the sensitivity median defined in [55].

Associated production of CP odd Higgs with heavy quarks
As a result of electroweak symmetry breaking and the Stueckelberg mass growth, the Higgs sector of the U (1) X -extended MSSM has six degrees of freedom corresponding to three CP even Higgs, h, H and ρ and one CP odd Higgs A and two charged Higgs H ± . In this section we discuss the production and decay of the CP odd Higgs A. Since the weak scale of supersymmetry is high lying in the few TeV region we are in the so-called decoupling limit where the light CP even Higgs h is SM-like and H, A and H ± (charged Higgs) have comparable masses and much greater than h. In this case, A exhibits no tree-level couplings to the gauge bosons and couplings to down (up) type fermions that are (inversely) proportional to tan β. For high tan β values, tan β 10, the H, A Yukawa couplings to bottom quarks and tau leptons are strongly enhanced, while those to top quarks are strongly suppressed. In this region, the b-quark will play an important role as its coupling to the CP odd Higgs is enhanced. For this reason, we examine the associated production of A with bottom-anti bottom quarks, bbA.
There are two approaches to calculating the production cross-section of A in association with bb. The first considers the b-quark to be heavy and appearing only in the final state as shown in Fig. 2 where the leading order (LO) partonic processes are gg → bbA, qq → bbA.
(39) Figure 2: A sample of the tree (top two) and one-loop (bottom three) Feynman diagrams for bbA production at the LHC in the four-flavor scheme.
This approach constitutes the four-flavor scheme (4FS) in which the mass of the b-quark is considered part of the hard scale of the process. Apart from its dependence on the CP odd Higgs mass and tan β, the cross-section is sensitive to the b-quark mass which is taken to be the running mass at the appropriate renormalization and factorization scales. The LO 2 → 3 diagrams in the 4FS begin at O(α 2 S ) while the next-to-leading order (NLO) diagrams (bottom three diagrams of Fig. 2) contain bottom and top quarks circulating in the loops. Here the LO cross-section is proportional to α 2 S y 2 b while at NLO, the cross-section is proportional to α 3 S y 2 b and α 3 S y b y t , where y b , y t are the bottom and top Yukawa couplings and where the y b y t term corresponds to interference between the gluon-fusion (with a top quark in the loop) and bbA processes. As mentioned before, the CP odd Higgs coupling to the top quark is suppressed and thus the diagrams involving top quarks do not contribute significantly to the cross-section. Following the prescription of [59], the hard scale of the process is chosen at the renormalization and factorization scales such that µ R = µ F = (m A + 2m b )/4 with m A the CP odd Higgs mass and m b being the b-quark pole mass while the running b-quark mass is m b (µ F ) (see Table 4). The 4FS NLO cross-section at fixed order in α S is calculated with MadGraph5 aMC@NLO-2.6.3 using FeynRules [60] UFO files [61,62] for the Type-II two Higgs doublet model (2HDM). The choice of the latter is justified due to the fact that SUSY-QCD effects for our benchmarks are very minimal since the squarks and gluinos and heavy. The cross-sections at 14 TeV and 27 TeV are displayed in Table 4 along with uncertainties arising from scale variations.   Table 1. The running b-quark mass, in GeV, is also shown evaluated at the factorization and normalization scales, µ F = µ R (in GeV).
At any order in perturbation, the 4FS cross-section involves terms ∼ α n S log n (µ F /m b ) resulting from collinear splitting of gluons to bb pairs. This term is kept under control as long as µ F ∼ m b , however, this is not the case for our benchmarks especially for larger masses of the CP odd Higgs where such a term spoils perturbative convergence. The way to resolve this issue is by absorbing those terms to all orders in α S . The resummation of those potentially large logarithms is done via the DGLAP evolution of b-quark PDFs which constitutes the second approach to calculating cross-sections which is the five-flavor scheme (5FS). In this scheme, b-quarks are massless and considered as partons, so they do not appear in the final states at the partonic level. Hence the LO process (zeroth order in α S ) in the 5FS for CP odd Higgs production is bb → A.
At the parton level, the 4FS approach has advantage over that of 5FS since realistic btagging can be done with the former while the latter does not possess this property due to less rich final states. However, the 5FS parton level events are matched to parton showers which add b jets allowing proper b-tagging at the analysis level. This is of course pertinent to LO calculations while at higher orders in QCD, the 5FS start to exhibit richer final states with the appearance of b-quarks. The 5FS bbA production cross-section is known at next-to-NLO (NNLO) and we use SusHi-1.7.0 [63] to determine those cross-sections at 14 TeV and 27 TeV. The renormalization and factorization scales are µ R = m A and µ F = m A /4, respectively, which have been shown to be the suitable choices. Scale uncertainties are determined by varying µ R and µ F such that µ R , 4µ F ∈ {m A /2, m A , 2m A } with 1/2 ≤ 4µ F /µ R < 2. Although the b-quark is massless, the bottom Yukawa coupling is non-zero and renormalized in the MS scheme. The LO cross-section in 5FS is proportional to y 2 b while N k LO is proportional to y 2 b α k S with y b y t terms vanishing order-by-order in perturbative QCD. In calculating the cross-sections for both 4FS and 5FS cases we have used PDF4LHC15 nlo mc and PDF4LHC15 nnlo mc [64] PDFs, respectively. In order to combine both estimates of the cross-section, we use the Santander matching criterion [65] such that The matched cross-section of the inclusive process lies between the 4FS and 5FS values but closer to the 5FS value owing to the weight α which depends on the CP odd Higgs mass. The uncertainties are combined as such, In Table 4 we give the NLO 4FS, NNLO 5FS and matched cross-sections at 14 TeV and 27 TeV for the ten benchmarks of Table 1 along with µ R , µ F and the running b-quark mass in the 4FS case. Notice the dramatic increase in cross-section in going from 14 TeV to 27 TeV due to the production of strongly interacting particles along with A. The cross-sections have been checked with publically available results [59] by a proper scaling of the bottom Yukawa coupling. In the MSSM, the tree-level Higgs Yukawa coupling to bottom quarks is given by for tan β 1, where v is the SM VEV. Besides QCD corrections, this Yukawa coupling receives SUSY-QCD corrections given, at one-loop level, by [66] where mb 1,2 and mg are the sbottom and the gluino masses. Taking this correction into consideration, one then needs to scale the SM bbh cross-section by the square of in order to obtain the MSSM cross-section. However, ∆ b is negligible for our benchmarks due to heavy gluinos and sbottoms. In this case, the scaling only requires multiplying by tan 2 β to which we have found reasonable agreement with our results.
Due to its enhanced coupling to bottom quarks, the CP odd Higgs preferentially decays to bb pair while the second largest branching ratio is to τ + τ − pair as shown in Table 5. In the MSSM, the branching ratio to Z h is quite small and is not considered as a significant channel for discovery, at least at 14 TeV.  , the results exclude tan β > 5 for m A = 250 GeV and tan β > 51 for m A = 1500 GeV. Due to the absence of light neutralinos in the spectrum, the hMSSM 2 provides more stringent constraints because of higher A → τ τ branching ratio. Thus here tan β > 1 for m A = 250 GeV and tan β > 42 for m A = 1500 GeV are excluded. CMS uses the same benchmarks, cross-sections and branching ratios of the CP odd Higgs and arrives at similar exclusion limits where m A 250 GeV is excluded for tan β > 6 and the exclusion contour reaches 1600 GeV for tan β = 60. Projections for HL-LHC studies regarding the mass reach for the CP odd Higgs in case no excess is found is available [73]. The benchmarks of Table 1 are not yet excluded by experiment and lie within the contour set for HL-LHC. It is worth stressing the fact that such interpretations of experimental results as mentioned above are carried out within models that are very different from the one we consider here.
The signal we investigate consists of a CP odd Higgs decaying to two hadronic taus and produced alongside two b-quarks which can be tagged. Even in the 5FS, b-flavored jets can appear at the parton shower level and so b-tagging is viable here too. In order to account for misidentified b-tagged jets, we require that our final states contain at least one b-tagged jet and two tau-tagged (τ h ) jets such that p T (b) > 20 GeV, |η(b)| < 2.5 and p T (τ h ) > 15 GeV.
The standard model backgrounds relevant to the final states considered here are tt, t+jets, t+W/Z, QCD multijet, diboson and W, Z/γ * +jets. The signal and SM backgrounds are simulated at LO with MadGraph5-2.6.3 interfaced with LHAPDF [74] and using the NNPDF30LO PDF set. The cross-sections of the SM backgrounds are then normalized to their NLO values while those of the signal are scaled to their matched values in Table 4. The parton level events are passed to PYTHIA8 [75] for showering and hadronization. A five-flavor MLM matching [76] is performed on the backgrounds to avoid double counting of jets at the shower level. Jets are clustered with FASTJET [77] using the anti-k t algorithm [78] with jet radius 0.4. Detector simulation and event reconstruction is handled by DELPHES-3.4.2 [79] using the new cards for HL-LHC and HE-LHC generic detectors. The resulting files are read and analyzed with ROOT-6.16 [80].
Due to the smallness of the signal cross-section in comparison to the SM backgrounds (especially following the selection criteria), we use Boosted Decision Trees (BDT) to separate the signal from the background. The type of BDT used here is known as "Adaptive BDT" (AdaBoost). Before giving a brief overview of BDTs, we list the kinematic variables used to help in discriminating the signal from the background: 1. The total transverse mass of the di-tau system is given by [81] where This variable has the best separating power especially for heavier CP odd Higgs mass.
2. The hadronic di-tau invariant mass, m τ h τ h , has the same effect as m tot T and in addition works well for low mass signals.
4. The number of charged tracks associated with the leading tau, N τ tracks . Due to its one and three-prong decays, a tau can be identified by the tracks' charge multiplicities.
5. Due to the presence of b-tagged jets, we use the number of such jets, N b jet , as a discriminating variable.
6. Due to the rich jetty final states, we define the variable ln(p jet T ) as where p jet 1 T is the p T of the leading jet.
7. The di-jet transverse mass m di−jet T of the leading and sub-leading jets. It is a good discriminant against QCD multijet which tends to have a large value of this variable. If two jets cannot be found in an event, this variable is set to zero.

The effective mass defined as
where H T is the sum of the hadronic p T 's in an event, p T (τ h1 ) and p T (τ h2 ) are the transverse momenta of the leading and sub-leading hadronic taus.
BDTs employ a multivariate analysis technique to a classification problem such as the one at hand. The aim here is to classify a certain set of events as belonging to the signal or the background by using a number of discriminating variables (1−8 listed above) to make the decision. The signal (S) and background (B) samples undergo two phases: the training and testing phases. In the first phase, the BDTs are trained on those samples using the available list of kinematic variables. The algorithm sorts those variables in a descending order of separation power and chooses the variable that has the highest separation power to start the "root node". A cut is applied on this variable and events are split into left or right nodes depending whether they were classified as signal or background. Afterwards, another variable (or sometimes the initiating variable is kept) is chosen with a cut value which further splits down events into signal and background. The tree continues growing until a stopping criterion, such as the tree depth, is reached. The end layer of the tree contains the leaves which host the events classified as signal and given a value +1 or background and given a value −1. During training, some signal events may be misclassified as background and vice-versa. Those events will be given a weight factor and then enter in the second iteration of the training phase when the next tree starts forming. Those events are now given more attention thanks to the weight factor they carry. The training stops when the entire number of trees in the forest have been utilized. The number of trees and their depths are specified by the user in such a way as to maximize the separation between the signal and the background where a larger depth generally produces a better separation.
The second phase is the testing phase where the algorithm applies what it has learnt on a statistically independent set of samples and outputs a new discriminating variable called the "BDT score" or "BDT response". An agreement between the performances of the training and testing phases is a sign of no overtraining occurring in the analysis. Such a situation can arise if one chooses too large of a tree depth while not having enough statistics in the samples. We have made sure that no overtraining of the samples has occurred throughout. The BDT implementation is carried out using ROOT's own TMVA (Toolkit for Multivariate Analysis) framework [82]. Depending on the samples, we set the number of trees to be in the range 120 to 200, the depth to 3 and the AdaBoost learning rate to 0.5. Many combinations of those parameters have been tried and the one which gave the best result was considered. BDTs are very useful in classification problems where conventional linear cuts fail. To show that, we display in Fig. 3 distributions normalized to unity of four kinematic variables at √ s = 27 TeV for benchmark (a) of Table 1. The purpose of such distributions is to help design event selection cuts which would allow better background rejection based on the shape of the distribution. One can clearly notice, for distributions in ln(p jet T ), N τ tracks and ∆φ(τ h1 , τ h2 ), a conventional linear cut does not do the job as it would lead to a poor signal to background ratio. This is where BDTs become powerful since they employ non-linear cuts in the multidimensional space of variables (thus the name multivariate analysis). On the other hand, a linear cut on m tot T such that m tot T > 150 GeV is reasonable but not sufficient to obtain a good signal to background ratio. The BDT algorithm will run through distributions of such sort for the eight variables presented above and design unconventional cuts to obtain the best discrimination between the signal and the background. After training and testing of the BDTs, we set a cut on the BDT score variable which would give us the minimum integrated luminosity for S √ S+B at the 5σ level discovery. In general, this cut value is not common across all points since each point is trained and tested separately along with the SM backgrounds and so the distribution in BDT score differs from one point to another. We present in Fig. 4 the computed integrated luminosities, L, as a function of the cut on the BDT score for both 14 TeV (left panel) and 27 TeV (right panel) machines. For 14 TeV, one can see that a drop in L occurs for BDT score > 0.3 while at 27 TeV the same is seen for BDT score > 0.2. Using the results from Fig. 4, we tabulate the lowest integrated luminosity for a 5σ discovery at HL-LHC and HE-LHC in Table 6 for all our benchmarks. It is seen that half the benchmarks are discoverable at HL-LHC. Thus benchmark (d) is discoverable with an L as low as 866 fb −1 while the benchmark (f) requires L close to the optimal integrated luminosity of 3000 fb −1 . However, all the benchmarks are discoverable at HE-LHC with some requiring an integrated luminosity smaller than 100 fb −1 such as point (d) with L = 50 fb −1 for discovery. Point (j) requires the largest amount of data at ∼ 2600 fb −1 which, however, is still much lower than the optimal integrated luminosity of 15 ab −1 expected at HE-LHC. ... 2599 Table 6: Comparison between the estimated integrated luminosity (L) in fb −1 for a 5σ discovery at 14 TeV (middle column) and 27 TeV (right column) for the CP odd Higgs following the selection criteria and BDT cut. An entry with an ellipsis means that the evaluated L is much greater than 3000 fb −1 .
We show in Fig. 5  It is worth mentioning that in using BDTs, it is seen that better separation between signal and background occurs for points with higher CP odd Higgs masses due to more energetic final states. However, a better outcome, i.e. a smaller integrated luminosity, is not always seen in those cases since the price to be paid for larger masses is a falling cross-section which results in much higher integrated luminosity for discovery. Here we investigate the possibility if at higher masses the cross-section, σ × BR(A → τ τ ), can be maintained at a larger than usual value. This can be achieved for a higher branching ratio and larger tan β. Benchmark (b) has the largest cross-section amongst all the points but requires the second lowest integrated luminosity for discovery with point (d) requiring the least. Now point (d) has a CP odd Higgs mass 100 GeV heavier than point (b) but has a higher tan β and branching ratio which makes up for the mass increase and keeps the cross-section from falling too rapidly. Because of the effect of tan β, the branching ratio and more energetic final states, the separation between signal and background for point (d) is more pronounced leading to the lowest integrated luminosity for discovery.
Given that the HE-LHC is expected to collect data at the rate of 820 fb −1 per year [83] versus the rate at which HL-LHC will collect data, the projected runtime for points (a)−(d) for discovery at HL-LHC is ∼ 3 to ∼ 4 years while point (f) requires ∼ 8 years. The runtime is significantly decreased for HE-LHC where most of the points require ∼ 1 to ∼ 6 months, while point (d) ∼ 22 days and point (j) ∼ 3 years.
Before concluding we give an overview of the uncertainties one might expect and their impact on the estimated integrated luminosities at HL-LHC and HE-LHC. One of the main challenges to overcome in experiment while analyzing data are the systematic uncertainties. One would expect such uncertainties to decrease when HL-LHC starts operation due to an increased data set. It is also reasonable to assume that improvements on this front is expected by experimentalists working on ATLAS and CMS detectors. Following experts' opinion in the HL-LHC and HE-LHC Working Groups at CERN [84,85], much of the systematic uncertainties are expected to drop by a factor of 2. Those uncertainties are known as "YR18 systematic uncertainties". In the A → τ h τ h channel, systematic uncertainties due to the estimation of QCD jets to τ h fake background are dominant especially in the low CP odd Higgs mass region. For higher masses, the leading uncertainty is from the reconstruction and identification of high transverse momentum τ h . In this analysis we assume that the systematic uncertainties in the background and signal are 20% and 10%, respectively. We give higher (lower) combined uncertainty for points with smaller (larger) Higgs mass. The integrated luminosity for a 5σ discovery is re-estimated after including the uncertainties using the signal significance [86] σ = 2 (S + B) ln (S + B)(B + ∆ 2 C ) where ∆ C is the combined uncertainty in signal and background, ∆ 2 C = ∆ 2 S + ∆ 2 B . We show in Fig. 6 the estimated integrated luminosities before and after including the uncertainties for both HL-LHC and HE-LHC. In the left panel, the five benchmarks discoverable at both machines are shown along with the "YR18 uncertainties" where at HL-LHC, the integrated luminosities have increased by ∼ 1.5 to ∼ 2.5 times (in blue) compared to when no systematic uncertainties are present (in orange). At the HE-LHC the increase is by ∼ 1.5 to ∼ 4 times (in red) compared to the case with no systematics (in yellow). The right panel shows the points that are discoverable only at HE-LHC along with the integrated luminosities before (in orange) and after (in blue) including uncertainties.  Table 1 that are discoverable at both HL-LHC and HE-LHC with and without the "YR18 uncertainties". Right panel: the remaining five benchmarks of Table 1 that are discoverable only at HE-LHC with and without the "YR18 uncertainties".

Conclusions
The large size of weak scale supersymmetry implied by the Higgs boson mass at ∼ 125 GeV has a direct implication for the discovery of supersymmetric dark matter. Thus typically in high scale SUGRA models with a large size of weak scale supersymmetry often the LSP neutralino is mostly a bino making an efficient annihilation of bino dark matter in the early universe difficult and leading to its overabundance inconsistent with observation. Often coannihilation is utilized in this case where a low-lying next-to-LSP and the LSP together act to reduce the relic density within the observed limit. However, another branch of radiative breaking of the electroweak symmetry exists within the SUGRA model where a large size of weak scale supersymmetry can coexist with a small Higgs mixing parameter µ (of size the electroweak scale). Such a situation can lead to a higgsino-like LSP. As mentioned in the introduction, models of this type are severely constrained by simultaneous satisfaction of dark matter relic density and by the spin-independent proton-DM scattering cross section limits in direct detection experiments. However, such models can be made viable if the dark matter is multi-component. Thus in this work we considered an extension of the MSSM/SUGRA gauge group with a U (1) X of the hidden sector where the U (1) X and U (1) Y have kinetic and Stueckelberg mass mixings. Further, we assume that the hidden sector is populated with chiral matter leading to a Dirac fermion which acts as the second component of dark matter and makes up the dark matter deficit in the relic density. One implication of the model is the existence of a relatively light CP odd Higgs boson A, as well as relatively light H and H ± , which have masses of the electroweak size. We investigate a set of benchmarks for the extended model and show that the CP odd Higgs boson in models of this type is observable in the next generation collider experiments. Specifically we investigate the discovery potential of a CP odd Higgs in the τ h τ h final state at the HL-LHC and HE-LHC. It is seen that a CP odd Higgs with a mass up to 450 GeV for tan β ≤ 12 may be discoverable at HL-LHC. The discovery reaches 750 GeV at the HE-LHC with an integrated luminosity of ∼ 2600 fb −1 which is just a fraction of the optimal luminosity of 15 ab −1 that HE-LHC can deliver. With the optimal luminosity the mass reach of HE-LHC for the CP odd Higgs mass will certainly extend far above 750 GeV. It is also shown that a significant part of the parameter space of the extended model can be probed in the next generation direct detection experiments such as XENONnT and LUX-ZEPLIN.