1 Introduction

The Standard Model (SM) of fundamental interactions has been extensively tested in recent decades and the search for its last missing piece – the SM Higgs particle – ended in 2012 with the discovery of a scalar boson with a mass of approximately 125 GeV by ATLAS and CMS experiments at the CERN Large Hadron Collider (LHC) [1, 2]. Since then, further effort has been spared to study Higgs boson dynamics at the LHC. Although the properties of the observed scalar are in agreement with those of the SM Higgs boson, it is still possible that it is just one member of an extended (pseudo)scalar sector.

There are various reasons why it is generally believed that the SM of particle physics is incomplete. One of the issues that needs to be addressed is the absence of a Dark Matter (DM) candidate in the SM. Cosmological observations imply that about 85% of matter in the Universe is cold (i.e., non-relativistic at the onset of galaxy formation), non-baryonic, neutral and weakly interacting [3]: such a state does not exist in the SM. Various candidates have been proposed so far, the best studied being a Weakly Interacting Massive Particle (WIMP) [4,5,6]. The mass of this hypothetical particle can vary between a few GeV and a few TeV, however, its exact nature is still unknown.

A particle with such characteristics can come from an extended scalar sector with a discrete symmetry. A well-known example is the Inert Doublet Model (IDM), a 2-Higgs Doublet Model (2HDM) with an unbroken discrete \(Z_2\) symmetry [7]. The model involves 1 inert doublet, which is \(Z_2\)-odd, does not develop a Vacuum Expectation Value (VEV) and – by construction – does not couple to fermions, plus 1 active \(Z_2\)-even Higgs doublet, which has a non-zero VEV and couples to fermions in the same way as the SM Higgs doublet. Therefore we shall also refer to the IDM as the I(1+1)HDM to explicitly show the number of inert (I) and active Higgs (H) doublets. An important feature of this model is that, due to the unbroken \(Z_2\) symmetry, the lightest neutral \(Z_2\)-odd particle, coming from the inert doublet, is stable and a suitable DM candidate.

The I(1+1)HDM, despite being severely constrained by data, remains a viable model for a scalar DM candidate (see the latest analyses, e.g., in [8,9,10,11]). This model, by construction, can not contain CP-violation: due to the presence of an exact \(Z_2\) symmetry, all parameters in the potential are real. In fact, accommodating CP-violation in multi-inert models requires at least three scalar SU(2) doublets, leading to a 3-Higgs Doublet Model (3HDM). Here, one can have two possibilities.

  • I(1+2)HDM: a 3HDM with 1 inert doublet plus 2 active Higgs doublets,

  • I(2+1)HDM: a 3HDM with 2 inert doublets plus 1 active Higgs doublet.

In the I(1+2)HDM, the inert sector is identical to that of the I(1+1)HDM and CP-violation is introduced in the extended active sector [12, 13]. Therefore, the amount of CP-violation is restricted by SM Higgs data, as the Higgs particle observed at the LHC is very SM-like, and by contributions to the Electric Dipole Moments (EDMs) of electron and neutron.

In the I(2+1)HDM, in contrast, the active sector is by construction SM-like, with tree-level interactions identical to those of the SM Higgs, with the exception of possible Higgs decays to new states provided they are sufficiently light.Footnote 1 Here, the inert sector is extended and now contains six new particles, four neutral and two charged ones, i.e., twice as many inert particles as in the I(1+1)HDM. As a result, even without introducing CP-violation, the I(2+1)HDM provides new coannihilation channels for the DM candidate and revives regions of parameter space that are excluded in the I(1+1)HDM [14, 16]. With the introduction of CP-violation in the inert sector, the neutral inert particles will have a mixed CP quantum number. Note that the inert sector is protected by a conserved \(Z_2\) symmetry from coupling to the SM particles, therefore, the amount of CP-violation introduced here is not constrained by EDM data. The DM candidate, in this scenario, is the lightest state amongst the CP-mixed inert states which enlivens yet another region of viable DM mass range, with respect to both I(1+1)HDM and CP-conserving I(2+1)HDM [17].

In this paper, we study electron-positron collider signatures of a CP-violating I(2+1)HDM via the process \(e^+e^- \rightarrow Z^* \rightarrow S_iS_j\) (\(i,j=1, \ldots , 4\)), which has six possible final states, \(S_1S_{2,3,4}\), \(S_2S_{3,4}\), \(S_3S_4\) in the CP-violating case, in comparison to four possible final states, \(H_1 A_{1,2}, H_2A_{1,2}\) in the CP-conserving case, wherein \(H_{1,2}\) and \(A_{1,2}\) have opposite CP-parity.Footnote 2 Hence, a simple collider energy scan combined with a trivial counting experiment in the detectors revealing six thresholds rather than four will be a clear evidence of CP-violation, whether or not such \(S_i\) states will have been previously discovered.Footnote 3 This signal by itself does not provide a conclusive evidence for CP-violation, as the observable is a CP-even quantity. However, provided we have observed other processes hinting at a 3HDM, as in our previous and upcoming publications [14,15,16,17,18, 46], this signal will eventually help to distinguish between a CP-violating and a CP-conserving 3HDM by considering the number of observable Z boson decay thresholds in the pair production of neutral scalar states. In order to study this phenomenology, we provide several Benchmark Points (BPs), in agreement with all experimental and theoretical bounds, for which we show that the cross section of the \(e^+e^- \rightarrow Z^* \rightarrow S_iS_j\) process could be as large as a few picobarns at \(\sqrt{s}\) values accessible by future \(e^+e^-\) colliders. The proximity (or otherwise) of these thresholds would serve as characteristic signatures of different BPs with different DM properties.

The layout of the remainder of this paper is as follows. In Sect. 2, we present the details of the scalar potential and the theoretical and experimental limits on its parameters. In Sect. 3, we construct and justify our BPs. In Sect. 4, we show the production cross sections and decay thresholds in our BPs. In Sect. 5, we conclude and present the outlook for our future studies.

2 The scalar sector of the I(2+1)HDM

A 3HDM potential symmetric under a group G of phase rotations can be divided into two parts: a phase invariant part, \(V_0\), and a collection of extra terms ensuring the symmetry group G, \(V_G\) [19, 20]. Here, we consider a \(Z_2\)-symmetry, under which the three Higgs doublets \(\phi _{1,2,3}\) transform, respectively, as:

$$\begin{aligned} g_{Z_2}= \mathrm {diag}\left( -1, -1, 1 \right) . \end{aligned}$$
(1)

The resulting potential is of the following formFootnote 4:

$$\begin{aligned} V_{3HDM}= & {} V_0+V_{Z_2}, \nonumber \\ V_0= & {} - \mu ^2_{1} (\phi _1^\dagger \phi _1) -\mu ^2_2 (\phi _2^\dagger \phi _2) - \mu ^2_3(\phi _3^\dagger \phi _3) \nonumber \\&+\, \lambda _{11} (\phi _1^\dagger \phi _1)^2+ \lambda _{22} (\phi _2^\dagger \phi _2)^2 + \lambda _{33} (\phi _3^\dagger \phi _3)^2 \nonumber \\&+\, \lambda _{12} (\phi _1^\dagger \phi _1)(\phi _2^\dagger \phi _2) + \lambda _{23} (\phi _2^\dagger \phi _2)(\phi _3^\dagger \phi _3)\nonumber \\&+ \,\lambda _{31} (\phi _3^\dagger \phi _3)(\phi _1^\dagger \phi _1) \nonumber \\&+\, \lambda '_{12} (\phi _1^\dagger \phi _2)(\phi _2^\dagger \phi _1) + \lambda '_{23} (\phi _2^\dagger \phi _3)(\phi _3^\dagger \phi _2)\nonumber \\&+\, \lambda '_{31} (\phi _3^\dagger \phi _1)(\phi _1^\dagger \phi _3), \nonumber \\ V_{Z_2}= & {} -\mu ^2_{12}(\phi _1^\dagger \phi _2)+ \lambda _{1}(\phi _1^\dagger \phi _2)^2 + \lambda _2(\phi _2^\dagger \phi _3)^2 \nonumber \\&+\, \lambda _3(\phi _3^\dagger \phi _1)^2 + h.c. \end{aligned}$$
(2)

The parameters of \(V_0\) are by construction real. We allow for the parameters of \(V_{Z_2}\) to be complex, hence introducing explicit CP-violation in the model.

The doublets are defined as

$$\begin{aligned} \phi _1= \left( \begin{array}{c}{ H^+_1 } \\ \frac{H_1+iA_1}{\sqrt{2}} \end{array}\right) ,\quad \phi _2= \left( \begin{array}{c}{ H^+_2 } \\ \frac{H_2+iA_2}{\sqrt{2}} \end{array}\right) , \quad \phi _3= \left( \begin{array}{c}{ G^+ } \\ \frac{v+h+iG^0}{\sqrt{2}} \end{array}\right) , \end{aligned}$$
(3)

where \(\phi _1\) and \(\phi _2\) are the two inert doublets, \(\langle \phi _1 \rangle = \langle \phi _2 \rangle =0\), while \(\phi _3\) is the one active doublet, \(\langle \phi _3 \rangle =v/\)\( \sqrt{2} \) \( \ne 0\), and plays the role of the SM Higgs doublet, with h being the SM Higgs boson and \(G^\pm ,~ G^0\) the would-be Goldstone bosons.

We assign \(Z_2\) charges to each doublet according to the \(Z_2\) generator in Eq. (1): odd-\(Z_2\) charge to the inert doublets, \(\phi _1\) and \(\phi _2\), and even-\(Z_2\) charge to the active doublet, \(\phi _3\). It is clear that the symmetry of the potential is respected by the vacuum alignment (0, 0, v/\( \sqrt{2} \)). To make sure that the entire Lagrangian and not only the scalar potential is \(Z_2\) symmetric, we assign an even \(Z_2\) parity to all SM particles, identical to the active doublet \(\phi _3\). With this parity assignment Flavour Changing Neutral Currents (FCNCs) are avoided as the extra doublets are forbidden to couple to fermions and, as dictated by the \(Z_2\) symmetry, \(\phi _3\) is the only doublet that couples to the fermions though Yukawa interactions identical to those in the SM Yukawa Lagrangian:

$$\begin{aligned} \mathcal {L}_\mathrm{Yukawa}= & {} \Gamma ^u_{mn} \bar{q}_{m,L} \tilde{\phi }_3 u_{n,R} + \Gamma ^d_{mn} \bar{q}_{m,L} \phi _3 d_{n,R} \nonumber \\&+ \, \Gamma ^e_{mn} \bar{l}_{m,L} \phi _3 e_{n,R} + \Gamma ^{\nu }_{mn} \bar{l}_{m,L} \tilde{\phi }_3 {\nu }_{n,R} + h.c. \end{aligned}$$
(4)

Here, \(\Gamma ^{u,d,e,\nu }_{mn}\) are the dimensionless Yukawa couplings for the family indices m, n and \(u,d,e,\nu \) labels refer to the SM fermions in the usual notation.

Note that the scalar h contained in the doublet \(\phi _3\) in our model has the tree-level couplings of the SM Higgs boson. Thus CP-violation is only introduced in the inert sector which is forbidden from mixing with the active sector by the \(Z_2\) symmetry, so that the amount of CP-violation is not limited by EDMs. The lightest amongst the neutral fields from the inert doublets, which now have a mixed CP-charge, \(S_1, S_2, S_3, S_4\), is the DM candidate, indeed stable due to the unbroken \(Z_2\) symmetry. (We avoid regions of parameter space where one of the charged inert scalars is the lightest.)

2.1 The parameters of the potential

The parameters of the potential can be divided into the following categories.

  • The Higgs sector parameters

    \(\mu ^2_3, \lambda _{33}\) are Higgs field parameters, fixed by the Higgs mass. We use the value 125 GeV for the latter and from extremum conditions we have:

    $$\begin{aligned} m^2_h = 2\mu ^2_3 = 2\lambda _{33} v^2. \end{aligned}$$
    (5)

  • The dark sector parameters

    \(\lambda _1, \lambda _{11},\lambda _{22},\lambda _{12}, \lambda '_{12}\) are inert/dark sector parameters (inert scalars self-interactions) and in tree-level analysis they are only constrained through perturbative unitarity and positivity of V. Apart from that, they do not play any role in our analysis, as they do not influence tree-level DM and collider phenomenology. We therefore set them to a fixed value of 0.1.

  • The phenomenologically relevant parameters

    \(\mu ^2_{1},\mu ^2_{2},\mu ^2_{12}, \lambda _{31},\lambda _{23},\lambda '_{31},\lambda '_{23},\lambda _{2}, \lambda _{3}\) are related to masses of inert scalars and their couplings with the visible sector. These 9 parameters can in principle be determined by independent masses, mixing angles or couplings and the ranges that we allow for them in our numerical studies are

    $$\begin{aligned}&-10~ \text{ TeV }^2< \mu ^2_{1},\mu ^2_{2},\mu ^2_{12}< 10~ \text{ TeV }^2, \nonumber \\&\quad -0.5< \lambda _{31},\lambda _{23},\lambda '_{31},\lambda '_{23},\lambda _{2}, \lambda _{3} < 0.5. \end{aligned}$$
    (6)

    The only parameters here that can be complex are \(\mu ^2_{12}\), \(\lambda _2\) and \(\lambda _3\) for which we use the following notation

    $$\begin{aligned} \mu ^2_{12}= & {} \mathrm {Re }\mu ^2_{12} +i \mathrm {Im }\mu ^2_{12} = |\mu ^2_{12}| e^{i \theta _{12}},\nonumber \\ \lambda _2= & {} \mathrm {Re }\lambda _2 +i \mathrm {Im }\lambda _2 = |\lambda _2| e^{i \theta _2}, \nonumber \\ \lambda _3= & {} \mathrm {Re }\lambda _3 +i \mathrm {Im }\lambda _3 = |\lambda _3| e^{i \theta _3}. \end{aligned}$$
    (7)

    Note that the phase of \(\mu ^2_{12}\) is non-physical and can be rotated away with the following redefinition of doublets

    $$\begin{aligned}&\phi _1 \rightarrow \phi _1 e^{i \theta _{12}/2} ~ \quad \quad |\mu ^2_{12}| e^{i \theta _{12}} \rightarrow |\mu ^2_{12}| , \nonumber \\&\phi _2 \rightarrow \phi _2 e^{-i \theta _{12}/2} \Longrightarrow |\lambda _2| e^{i \theta _2} \rightarrow |\lambda _2| e^{i (\theta _2+\theta _{12})}, \nonumber \\&\phi _3 \rightarrow \phi _3 ~ \quad \quad \quad \quad |\lambda _3| e^{i \theta _3} \rightarrow |\lambda _3| e^{i (\theta _3+\theta _{12})}. \end{aligned}$$
    (8)

    We, therefore, set \(\theta _{12}\) to zero for simplicity.

2.1.1 The dark democracy limit

In our previous papers [14, 16, 17], we studied a simplified version of the I(2+1)HDM by imposing the following equalities

$$\begin{aligned} \mu ^2_1 =\mu ^2_2 , \quad \lambda _3=\lambda _2 , \quad \lambda _{31}=\lambda _{23} ,\quad \lambda '_{31}=\lambda '_{23} , \end{aligned}$$
(9)

which is sometimes referred to as the dark democracy limit. After imposing this limit, the model is still explicitly CP-violating when \((\lambda _{22}- \lambda _{11} ) [\lambda _1 ({\mu ^2_{12}}^*)^2-\lambda ^*_{1}(\mu ^2_{12})^2 ] \ne 0\) [21, 22]. Note that, after rotating away the phase of \(\mu _{12}^2\), the amount of CP violation is directly related to the dark sector through the parameters \(\lambda _{11,22}\) and a complex \(\lambda _1\).

2.1.2 The dark hierarchy limit

In this paper, we study the more general case of the darkhierarchy:

$$\begin{aligned}&\mu ^2_1 =n\mu ^2_2 , \quad \text {Re}\lambda _3=n \text {Re}\lambda _2 , \quad \text {Im}\lambda _3=n \text {Im}\lambda _2 ,\quad \nonumber \\&\lambda _{31}=n \lambda _{23} ,\lambda '_{31}=n \lambda '_{23} , \end{aligned}$$
(10)

where we introduce the dark hierarchy parameter n, which can change between \(0 \le n \le 1\). Boundary values reduce the model to the well-known I(1+1)HDM for \(n=0\) and to the dark democracy case for \(n=1\). The case of \(n>1\) corresponds to a redefinition of states and does not lead to any different phenomenology.

After imposing the dark hierarchy limit, the only two relevant complex parameters, \(\lambda _2\) and \(\lambda _3\), are related through \(|\lambda _3 | = n |\lambda _2|\) and \( \theta _3 = \theta _2\). The angle \(\theta _2\) is therefore the only relevant CP-violating phase and is referred to as \(\theta _\mathrm{CPV}\) throughout the paper.

Note that n is a discrete parameter which we take as an input in our analysis. In reality, this quantity is subject to RGE effects meaning that our set-up may need readjustment in higher orders. Nonetheless, our low energy phenomenological analysis is always attainable by smooth changes of the n parameter.

2.2 Physical scalar states

The \(Z_2\)-conserving minimum of the potential sits at the point \((0,0,\frac{v}{\sqrt{2}})\) with \( v^2= \frac{\mu ^2_3}{\lambda _{33}}\). The resulting mass spectrum of the scalar particles is as follows.

2.2.1 The fields from the active doublet

The fields from the third doublet, \(G^0,G^\pm ,h\), which play the role of the SM Higgs doublet fields have squared masses of

$$\begin{aligned}&m^2_{G^0}= m^2_{G^\pm }=0, \nonumber \\&m^2_{h}= 2\mu _3^2 =2\lambda _{33} v^2. \end{aligned}$$
(11)

2.2.2 The charged inert fields

The two physical charged states, \(S_{1}^{\pm }\) and \(S_{2}^{\pm }\), from the inert doublets are the eigenstates of the matrix

$$\begin{aligned} \mathcal {M}_C= \left( \begin{array} {cc} -n \mu _{2 }^{2} + \frac{n}{2} \lambda _{23} v^{2} &{}\quad -\mu _{12}^{2} \\ - \mu _{12}^{2} &{}\quad -\mu _{2}^{2} + \frac{1}{2} \lambda _{23} v^{2} \end{array} \right) , \end{aligned}$$
(12)

with eigenvalues:

$$\begin{aligned} m^2_{S^\pm _{1,2}}= & {} \frac{1}{4} \biggl ((n+1)(-2 \mu _2^2+\lambda _{23}v^2)\; \nonumber \\&\mp \,\sqrt{16 (\mu _{12}^2)^2+(n-1)^2 \left( \lambda _{23} v^2-2 \mu _2^2\right) ^2}\biggr ). \end{aligned}$$
(13)

In terms of gauge states from Eq. (3) \(S^\pm _i\) are defined through:

$$\begin{aligned}&\left( \begin{array} {c} S_1^\pm \\ S_2^\pm \end{array} \right) = \left( \begin{array} {cc} \cos \alpha _c &{}\quad \sin \alpha _c \\ -\sin \alpha _c &{}\quad \cos \alpha _c \end{array} \right) \left( \begin{array} {cccc} H_1^\pm \\ H_2^\pm \end{array} \right) \quad \text{ with } \quad \tan 2\alpha _c \nonumber \\&\quad = \frac{2 \mu _{12}^2}{(n-1) (\mu _2^2 - \lambda _{23} v^2/2)}. \end{aligned}$$
(14)

We require \(\pi /2< \alpha _c < \pi \), so that \(m_{S_1^\pm } < m_{S_2^\pm }\).

2.2.3 The neutral inert fields

The neutral mass-squared matrix in the \((H_1,H_2,A_1,A_2)\) basis is

$$\begin{aligned} \mathcal {M}_N= \frac{1}{4}\left( \begin{array} {cccc} n \;\Lambda ^+_c &{}\quad -2\mu ^2_{12} &{}\quad -n \;\Lambda _s &{}\quad 0 \\ -2\mu ^2_{12} &{}\quad \Lambda ^+_c &{}\quad 0 &{}\quad \Lambda _s \\ -n \;\Lambda _s &{}\quad 0 &{}\quad n \;\Lambda ^-_c &{}\quad -2\mu ^2_{12} \\ 0 &{}\quad \Lambda _s &{}\quad -2\mu ^2_{12} &{}\quad \Lambda ^-_c \end{array} \right) , \end{aligned}$$
(15)

with

$$\begin{aligned} \Lambda _s= & {} 2\lambda _2 \sin \theta _\mathrm{CPV} v^2 \quad \text{ and } \quad \Lambda ^\pm _c =-2\mu ^2_2 + (\lambda _{23}+\lambda '_{23} \nonumber \\&\pm \, 2 \lambda _2 \cos \theta _\mathrm{CPV}) v^2. \end{aligned}$$
(16)

Note that, in the CP-conserving limit, \(\theta _\mathrm{CPV} = 0,\pi \) leads to \(\Lambda _s=0\) which reduces \(\mathcal {M}_N\) to a block diagonal matrix with no mixing between the states with opposite CP-parity, namely between \(H_{1,2}\) and \(A_{1,2}\).

We diagonalise the neutral mass-squared matrix numerically, \(\mathcal {M}_N^\mathrm{diag} = R^T \mathcal {M}_N R\), to derive our mass eigenstates, \(S_i\), in terms of the gauge eigenstates in Eq. (3),

$$\begin{aligned} \left( \begin{array} {cccc} S_1\\ S_2 \\ S_3 \\ S_4 \end{array} \right) = R_{ij} \left( \begin{array} {cccc} H_1\\ H_2\\ A_1 \\ A_2 \end{array} \right) . \end{aligned}$$
(17)

We adopt a notation where \(m_{S_1}< m_{S_2}<m_{S_3} <m_{S_4}\), hence choosing \(S_1\) as DM candidate. We use

$$\begin{aligned} |\mu ^2_{12}| , \;\lambda _{23}, \; \lambda '_{23}, \; \mu ^2_2, \; \lambda _2, \; \theta _\mathrm{CPV}, n \end{aligned}$$
(18)

as the set of input parameters to define our BPs in a forthcoming section.

2.3 Constraints on the parameters

In this section, we discuss the latest theoretical and experimental constrains that are applicable to our studies. The I(2+1)HDM is a model which is partially already within reach of current collider as well as DM experiments, and their results constrain parts of parameter space.

2.3.1 Theoretical constraints

In the “dark hierarchy” limit, theoretical requirements of boundedness of the potential and positive-definiteness of the Hessian are taken into account. All couplings fulfil perturbative unitarity limits, i.e., they take absolute values \(\lambda _i\le \,4\,\pi \), as noted in Eq. (6). For detailed formulas see [14, 16, 17].

2.3.2 Experimental constraints

  • Higgs decays and signal strengths

    The latest measurement of the SM-like Higgs boson’s width gives \(\Gamma _\text {tot}\,=3.2^{+2.8}_{-2.2}\) MeV, with 95% CL upper limit of \(\Gamma _\text {tot} \le 9.16\) MeV [23]. In our model, the total width of the SM-like Higgs boson can be modified through two mechanisms. If inert scalars are light, \(m_{S_i} <m_h/2\), we can expect a measurable contribution to Higgs invisible decays. This sets strong limits on the Higgs-inert couplings in the light mass region. Furthermore, the partial decay \(\Gamma (h\rightarrow \gamma \gamma )\) can be significantly changed through the two charged inert scalar contributions, as new physics corrections are formally of the same order as the SM process. In this work we use the combined ATLAS and CMS Run I limit for the signal strength \(h\rightarrow \gamma \gamma \), \(\mu _{\gamma \gamma } = 1.14^{+0.38}_{-0.36}\) [24].Footnote 5 By construction, the Higgs particle is SM-like and couplings to gluons, massive gauge bosons and fermions are equal to the SM values.

    The aim of this paper is to present a way of testing the model at future linear colliders. However, before either CLIC or ILC are built, the analysis of current and future runs at the LHC are going to provide stronger constraints on the parameter space. Therefore, it is also important to establish if the model is in agreement with projected LHC exclusion limits. The diphoton signal strengths \(\mu \) for benchmarks A, B, C (presented in details in following sections) are \(\mu _A = \mu _B = 0.937, \mu _C = 0.853\), respectively. They are all within the 3\(\sigma \) limit of current measurement. Note, that the signal strength is the same for two benchmark A and B, which correspond to a vastly different dark sector. This further stresses the importance of using linear colliders to fully explore the models with extended scalar sector. Assuming predictions from the HL-LHC [38], we expect \(\mu =1.00 \pm 0.04\) in the ggF channel, and \(\mu =1.00^{+ 0.10}_{-0.09}\) in the VBF channel, which means these benchmarks will still be in agreement with the observation.

    To be in agreement with current experimental constraints we have to minimize the contribution of new particles to Higgs decays. Two benchmarks have particles with masses below \(m_h/2\), benchmarks B and C. However, as \(\sum _{i,j} BR(h\rightarrow S_i S_j)_B = 0.07\%, \sum _{i,j} BR(h\rightarrow S_i S_j)_C = 0.5\%, i,j = 1,2,3,4\), we do not observe any significant changes in the total decay width of the Higgs with respect to the SM prediction. Furthermore, even taking into account the most promising estimation of the final constraint on Higgs invisible decays from HL-LHC and LeHC, \(\lesssim 4-7\%\) and \(\lesssim 2.25\%\), our chosen benchmarks will still be in agreement with the observation.

    In short, projected sensitivity limits from Higgs physics of future runs of the LHC are not going to provide any relevant constraints for the parameter range of the model considered in this paper.

  • Gauge bosons widths

    Similarly to the Higgs width, if new particles are sufficiently light, they could significantly change the total width of Electro-Weak (EW) gauge bosons. We control this by forbidding decays \(W^\pm \rightarrow S_i S_j^\pm \) and \(Z\rightarrow S_i S_j,S_i^+S_j^-\) through enforcing:

    $$\begin{aligned} m_{S_i}+m_{S^\pm _i}\,\ge \,m_W^\pm ,~~ \,m_{S_i}+m_{S_j}\,\ge \,m_Z,\,~~ 2\,m_{S_i^\pm }\,\ge \,m_Z. \end{aligned}$$
    (19)

  • EW Precision Observables (EWPOs)

    We require a \(2\sigma \), i.e., a \(95 \%\) Confidence Level (CL), agreement, parameterised through the EW oblique parameters S, T, U [27,28,29,30]. Just like in the 2HDM, it suffices here to have in the dark sector (near) degeneracy between each charged state and one or two of the neutral ones, condition which is satisfied by all our BPs.

  • Charged scalar mass and lifetime

    We take a conservative lower estimate on the masses of charged scalars [31] \(m_{S^\pm _i} > 70\) GeV (\(i=1,2\)). We also ensure that neither of these particles are quasi-stable by setting a limit for the charged scalar lifetime to be \(\tau \,\le \,10^{-7}\,\mathrm{s}\) [32].

  • Searches for new particles at colliders

    As in previous works of ours, we use LEP 2 searches for supersymmetric particles (chiefly, sneutrinos and sleptons) re-interpreted for the IDM in order to exclude the region of masses where the following condition are simultaneously satisfied [33] (\(i,j=2, \ldots , 4\)):

    $$\begin{aligned}&m_{S_i}\,\le \,100\,\mathrm {\;GeV},\,~~ m_{S_1}\,\le \,80\,\mathrm {\;GeV},\,\, ~~ \Delta \nonumber \\&m {(S_i,S_1)}\,\ge \,8\,\mathrm {\;GeV}, \end{aligned}$$
    (20)

    since this would lead to a visible di-jet or di-lepton signal.Footnote 6

    Benchmarks are in agreement with null-results for additional neutral scalar searches at the LHC, where we make use of HiggsBounds-5.4.0beta [39,40,41,42] and HiggsSignals-2.2.3beta [43]. This follows the analysis performed in [11] for the I(1+1)HDM, which is a model with similar signatures to one studied in this paper, especially for benchmark A.

    One of possible ways of testing the model would be using searches for multilepton final states with missing transverse energy. However, current strategy at the LHC uses a relatively large cut on missing transverse energy, which corresponds to a rather large mass splittings between scalars in the dark sector, and that reduces the production cross sections. Benchmarks with smaller mass splittings between scalars have large enough cross-section to be produced in abundance even at the current stage of the LHC, however they require smaller cuts on missing energy to be detected. With current analysis setup used at the LHC, the considered parameter range of the I(2+1)HDM model will remain invisible even for future HL-LHC or HE-LHC, especially for benchmarks of type B and C. We strongly encourage the LHC experimental collaborations to expand their search region in multilepton final states and missing \(E_T\) by allowing for lower cuts on missing energy, which would give us access to the considered parameter space of the model. Note that mono-photon searches at LEP may put more limiting constraints on light DM particles with light mediators [34, 35]. However, in our chosen benchmark points with \(m_{DM} \ge 50\) GeV and mediator mass of \(m_Z = 91\) GeV, direct detection is the more constraining probe of the model.

  • DM measurements

    We require agreement with relic density limits from the Planck experiment [3]:

    $$\begin{aligned} \Omega _c\,h^2\,=\,0.120\,\pm \,0.001, \end{aligned}$$
    (21)

    as well as with the latest XENON1T results for direct DM searches [36]. In the region of masses we are considering in this paper, indirect detection experiments (e.g., FermiLAT) do not place any additional constraints upon the parameter space.

    Note that, due to the absence of any \(S_iS_i Z\) vertex, the only loop-induced process that could potentially contribute to the DM-nucleon interaction is the rightmost diagram in Fig. 1. However, since direct detection experiments are performed at the zero-momentum limit, for this diagram to contribute to the DM-nucleon cross-section, the mass splitting between \(S_1\) and the next lightest inert scalar, \(S_{2,3,4}\), has to be of the keV order. In all our benchmark points, the \(S_i-S_j\) mass splitting is of above 1 GeV. Therefore, the only tree-level contribution to the DM-nucleon cross-section is from the \(S_1N\rightarrow S_1N\) process mediated by the Higgs boson.

    The one-loop contributions, representedFootnote 7 by the rightmost diagram in Fig. 1, are shown to have a relevant contribution in the I(1+1)HDM [37], which is CP conserving by construction. However, due to CP violation in our model, the strength of the gauge couplings is reduced. Therefore, the results of [37] are not directly applicable here. We postpone the detailed study of DM direct detection one-loop contributions to a future publication, since the focus of this paper is on the collider signatures of the model.

Fig. 1
figure 1

The only tree-level contribution to the direct detection cross-section is from the tree-level \(S_1N\rightarrow S_1N\) process mediated by the Higgs boson (left). With \(m_{S_{j}}-m_{S_1} \gg \) few keV (\(j=2,3,4\)), the \(S_1N\rightarrow S_jN\) processes do not contribute to the direct detection cross-section. The one-loop processes are represented by the rightmost diagram

Full size image

3 Selection of BPs

In our previous papers [14, 16, 17], we discussed DM phenomenology of the I(2+1)HDM in detail, where in addition to standard Higgs/gauge mediated annihilation channels of DM, there exist the possibility of coannihilation with heavier states, provided they are close in mass. This is a feature of models with extended dark sectors and contributes to changes in DM relic density. It is important to note that the relevance of this effect will depend not only on the DM mass and the mass splittings but also on the strength of the standard DM annihilation channels. For example, in some BPs presented in later sections, coannihilation, although possible, is responsible for less than 1% of the overall contributions to \(\Omega _\mathrm{DM}h^2\) because of a very strong DM annihilation into gauge bosons.

As detailed in [14, 16, 17], the generic expected behaviour of the I(2+1)HDM in different DM mass regions, is as follows.

  • Regions with \(m_\mathrm{DM} \le 45\) GeV are ruled out due to the Z gauge boson width constraints.

  • In the 45 GeV \( \le m_\mathrm{DM} \le 53\) GeV range, \(S_1\) mainly (co)annihilates with SM fermions,

    $$\begin{aligned} S_i S_j \rightarrow h / Z \rightarrow f \bar{f}. \end{aligned}$$
    (22)

    In this region, the \(h\rightarrow { invisible}\) channel is open and requires a very small Higgs-DM coupling to satisfy the experimental bounds.

  • In the 53 GeV \(\le m_\mathrm{DM} \le 75\) GeV range, (co)annihilations could also be mediated by the SM gauge bosons, \(V=Z, W^\pm \),

    $$\begin{aligned} S_i S_j \rightarrow V V^* \rightarrow V f \bar{f}, \qquad S_i S_j \rightarrow V^* V^* \rightarrow f\bar{f} f \bar{f}. \end{aligned}$$
    (23)

    The \(h\rightarrow invisible\) channel is closed, however, strong bounds from direct and indirect detection experiments require a very small Higgs-DM coupling.

  • In the 75 GeV \( \le m_\mathrm{DM} \lesssim 375\) GeV range, \(S_1\) (co)annihilation with gauge bosons,

    $$\begin{aligned} S_i S_j \rightarrow h \rightarrow V V, \qquad S_i S_j \rightarrow V V, \end{aligned}$$
    (24)

    is so strong that the model may fail to provide 100% of the observed DM relic density (so that a second DM component may need to be invoked, albeit in a wider framework than our I(2+1)HDM).

  • In the \(m_\mathrm{DM} \gtrsim 375\) GeV range, coannihilations with \(S_j^\pm \),

    $$\begin{aligned} S_i S_j^\pm \rightarrow W^\pm \rightarrow f f', \end{aligned}$$
    (25)

    interfere destructively with (co)annihilation to gauge bosons. As a result, the model provides 100% of the observed DM relic density.

Taking all theoretical and experimental bounds (listed in Sect. 2.2) into account, we have devised a few benchmark scenarios which show interesting phenomenology. In this paper, we do not consider the heavy mass region for the DM candidate (\(m_\mathrm{DM} \gtrsim m_Z\)), due to the fact that the \(e^+e^-\) production cross section of the heavy inert scalars drops significantly with an increasing \(m_\mathrm{DM}\) value (since \(m_\mathrm{DM}\equiv m_{S_1}<m_{S_2}<m_{S_3}<m_{S_4}\)). These points could be tested if \(\sqrt{s}\) is increased beyond the maximum value that we will consider, of 500 GeV. Also, as the heavy mass region corresponds to a semi-degenerate spectrum (in order to satisfy DM relic density bounds), we are not expecting to see there the interesting signatures and separation of thresholds that can be detected for the medium mass region (\(45\hbox { GeV }\le m_\mathrm{DM} \le m_Z\)), as it will be discussed in Sect. 4.

In such a medium mass region, the I(2+1)HDM provides three distinctive types of benchmark scenarios. To avoid all exclusion limits, we require a very small Higgs-DM coupling, \(g_{h\mathrm{DM}} \simeq 10^{-3}\).

  • Scenario A

    This is a case with large mass splittings, of order 50 GeV or so, between the DM candidate and all other inert particles:

    $$\begin{aligned} m_{S_1} \ll m_{S_2}, m_{S_3}, m_{S_4}, m_{S^\pm _1}, m_{S^\pm _2}. \end{aligned}$$
    (26)

    In this scenario no coannihilation channels are present and therefore \(S_1\) only annihilates through the Higgs boson to other SM particles.

    In the 45 GeV \( \le m_\mathrm{DM} < 53\) GeV range, the tiny \(g_{h\mathrm{DM}}\) does not provide an efficient annihilation of DM and is therefore forbidden by relic density observations. Within the mass range 53 GeV \( \le m_\mathrm{DM} \le 75 \) GeV, this scenario could easily accommodate points with a very small \(g_{h\mathrm{DM}}\) and avoid all exclusion limits.

  • Scenario B

    This is a case with a small mass splitting, of order 20% of \(m_\mathrm{DM}\), between the DM and one inert neutral particle:

    $$\begin{aligned} m_{S_1} \sim m_{S_2} \ll m_{S_3}, m_{S_4}, m_{S^\pm _1}, m_{S^\pm _2}. \end{aligned}$$
    (27)

    In this scenario the DM can coannihilate with its only particle close in mass, \(S_2\). This choice also leads to a relatively small mass splitting between \(S_3\) and \(S_4\), and effectively separating the neutral sector into two groups, with each generation accompanied by a charged scalar.

    In the 45 GeV \(\le m_\mathrm{DM} < 53\) GeV range, due to the existence of the coannihilation channel, DM is under-produced. In the mass range 53 GeV \(\le m_\mathrm{DM} \le 75\) GeV, where the destructive interference with coannihilation to gauge bosons comes into play, this scenario could accommodate points with very small \(g_{h\mathrm{DM}}\) and 100% DM contribution.

    If one relaxes the re-interpreted Supersymmetric limits, discussed in Sect. 2.2, and allows for larger mass splittings between \(S_1\) and \(S_2\), the strength of the \(S_1 S_2\) coannihilation channel could be reduced. As a result, this scenario can provide points where \(S_1\) contributes to 100% of DM relic density in the whole 45 GeV \( \le m_\mathrm{DM} \le 75\) GeV range with very small \(g_{h\mathrm{DM}}\).

  • Scenario C

    This is a case with all neutral particles close in mass:

    $$\begin{aligned} m_{S_1} \sim m_{S_2} \sim m_{S_3} \sim m_{S_4} \ll m_{S^\pm _1} \sim m_{S^\pm _2}. \end{aligned}$$
    (28)

    In this scenario the DM can coannihilate with all other neutral inert particles. Charged scalars are considerably heavier and do not participate in the coannihilation.

    Across the whole low and medium mass range, this scenario under-produces DM, due to the small mass splittings of the neutral inert particles which in turn strengthens the coannihilation channels. Contrary to the previous case, this situation cannot be resolved by relaxing the re-interpreted Supersymmetric limits and allowing for larger mass splittings. This is due to the large number of the coannihilation channels. As a result, with a very small \(g_{h\mathrm{DM}}\), this scenario will not be able to contribute to 100% of the observed relic density.

4 The \( e^{+} e^{-} \rightarrow Z^{*} \rightarrow S_{i} S_{j}\) cross section

We calculate the \( e^{+} e^{-} \rightarrow Z^{*} \rightarrow S_{i} S_{j}\) cross section at tree-level as [44, 45]

$$\begin{aligned} \sigma _{S_{i}S_{j}}&= \frac{ \pi \; \alpha ^{2} \; s \; g_{ZS_iS_j}^2}{24 \; (s-m_{Z}^2)^{2} \; g^2} \nonumber \\&\times \left( \frac{ 8 \; {\sin \theta _W}^4 - 4 {\sin \theta _W}^{2} + 1}{{\sin \theta _W}^{4} {\cos \theta _W}^{6}}\right) f^{3} (x,y), \end{aligned}$$
(29)

with \(x= m_{S_i}^{2}/s\), \(y=m_{S_j}^{2}/s\) and the function

$$\begin{aligned} f(x,y)= \sqrt{ 1 + x^{2} + y^{2} - 2x -2y -2xy}. \end{aligned}$$
(30)

The \(ZS_iS_j\) couplings are defined according to Eq. (31). Needless to say, cross sections with lighter final states will peak earlier at smaller \(\sqrt{s}\)) while those with larger \(g_{ZSiSj}\) coupling will be larger.

Following the discussion in Sect. 3, we have chosen three representative BPs from each possible scenario in the medium DM mass region. For all BPs, we aim to have at least one set of masses with \(m_{S_i} + m_{S_j} < 250\) GeV, which should lead to at least one channel being fully accessible at the first stage of a future \(e^+e^-\) collider, the so-called ‘Higgs factory’ run.

Below, for each BP, we list the input parameters, i.e., masses of particles and all relevant couplings, following the convention:

$$\begin{aligned} \mathcal {L}_\mathrm{gauge}&\supset g_{ZS_iS_j} Z_\mu (S_i \partial ^\mu S_j - S_j \partial ^\mu S_i), \end{aligned}$$
(31)
$$\begin{aligned} \mathcal {L}_\mathrm{scalar}&\supset \frac{v}{2}g_{S_i S_i h} h S_i^2+ v g_{S_i S_j h} h S_i S_j \nonumber \\&+\, v g_{S_i^\pm S_j^\mp h} h S_i^\pm S_j^\mp . \end{aligned}$$
(32)

4.1 Benchmark A

The input parameters for Benchmark A are defined as

$$\begin{aligned} \begin{array}{c} n= 0.6, \quad \lambda '_{23}= -0.16, \quad \lambda _{23}= 0.29, \quad \lambda _2= 0.067,\\ \theta _\mathrm{CPV}=15 \pi /16, \quad \mu _2^2= -13800~\mathrm{GeV}^2, \\ \mu _{12}^2= 5050~\mathrm{GeV}^2, \end{array} \end{aligned}$$
(33)

which lead to the following masses for the dark particles:

$$\begin{aligned} \begin{array}{c} m_{S_1}=72.331~\mathrm{GeV}, \quad m_{S_2}=103.313~\mathrm{GeV}, \\ m_{S_1^\pm }=106.235~\mathrm{GeV},\\ m_{S_3}=129.467~\mathrm{GeV}, \quad m_{S_4}=155.178~\mathrm{GeV}, \\ m_{S_2^\pm }=157.588~\mathrm{GeV}, \end{array} \end{aligned}$$
(34)

and the following gauge couplings:

$$\begin{aligned} \begin{array}{c} g_{ZS_1S_2}=0.366, \quad g_{ZS_1S_3}=0.0397, \quad g_{ZS_1S_4}=0.0401,\\ g_{ZS_2S_3}=0.04006, \quad g_{ZS_2S_4}=0.0397, \quad g_{ZS_3S_4}=0.366. \end{array} \end{aligned}$$
(35)
Fig. 2
figure 2

The \( e^{+} e^{-} \rightarrow Z^{*} \rightarrow S_{i} S_{j}\) cross section for BP A with masses in GeV

Full size image
Fig. 3
figure 3

The \( e^{+} e^{-} \rightarrow Z^{*} \rightarrow S_{i} S_{j}\) cross section for BP B with masses in GeV

Full size image

A characteristic signature of type A BPs, as shown in Fig. 2, is a pattern of very distinct thresholds that open up as \(\sqrt{s}\) increases, all easily resolvable thanks to the fine beam resolution available at future electron-positron machines. Here, the lightest and heaviest final states dominate over those with intermediate rest masses since the size of the cross section is dictated by the \(ZS_iS_j\) couplings.

4.1.1 Benchmark B

The input parameters for Benchmark B are defined as

$$\begin{aligned} \begin{array}{c} n=0.5, \quad \lambda '_{23}=-0.145, \quad \lambda _{23}=0.171, \\ \lambda _2= 0.013,\\ \theta _\mathrm{CPV}=7\pi /8, \quad \mu _2^2= -15900~\mathrm{GeV}^2, \\ \mu _{12}^2=7950~\mathrm{GeV}^2, \end{array} \end{aligned}$$
(36)

which lead to the following masses for the dark particles:

$$\begin{aligned} \begin{array}{c} m_{S_1}=55.441~\mathrm{GeV}, \quad m_{S_2}=63.219~\mathrm{GeV},\\ m_{S_1^\pm }=79.184~\mathrm{GeV},\\ m_{S_3}=144.377~\mathrm{GeV},\quad m_{S_4}=148.842~\mathrm{GeV},\\ m_{S_2^\pm }=159.203~\mathrm{GeV}, \end{array} \end{aligned}$$
(37)

and the following gauge couplings:

$$\begin{aligned} \begin{array}{c} g_{ZS_1S_2}=0.37, \quad g_{ZS_1S_3}=0.007, \quad g_{ZS_1S_4}=0.007,\\ g_{ZS_2S_3}=0.007, \quad g_{ZS_2S_4}=0.007, \quad g_{ZS_3S_4}=0.37. \end{array} \end{aligned}$$
(38)

A characteristic signature of type B BPs, as shown in Fig. 3, is two distinct thresholds, one (single) at low \(\sqrt{s}\) and another (single) at high \(\sqrt{s}\), plus several (very closely spaced) ones at mid \(\sqrt{s}\). Here too it is the lightest and heaviest final states that dominate over those with intermediate rest masses as dictated by the \(ZS_iS_j\) coupling strengths.

4.1.2 Benchmark C

The input parameters for Benchmark C are defined as

$$\begin{aligned} \begin{array}{c} n=0.8, \quad \lambda '_{23}=-0.295, \quad \lambda _{23}=0.294 \quad \lambda _2= 0.0009\\ \theta _\mathrm{CPV}=31 \pi /32, \quad \mu _2^2= -3400~\mathrm{GeV}^2, \quad \mu _{12}^2=250~\mathrm{GeV}^2, \end{array} \end{aligned}$$
(39)

which lead to the following masses for the dark particles:

$$\begin{aligned} \begin{array}{c} m_{S_1}=50.925~\mathrm{GeV},\quad m_{S_2}=51.793~\mathrm{GeV},\\ m_{S_1^\pm }=99.176~\mathrm{GeV},\\ m_{S_3}=58.555~\mathrm{GeV}, \quad m_{S_4}=59.459~\mathrm{GeV},\\ m_{S_2^\pm }=111.136~\mathrm{GeV}, \end{array} \end{aligned}$$
(40)

and the following gauge couplings:

$$\begin{aligned} \begin{array}{c} g_{ZS_1S_2}=0.37, \quad g_{ZS_1S_3}=0.0025, \quad g_{ZS_1S_4}=0.0028,\\ g_{ZS_2S_3}=0.0028, \quad g_{ZS_2S_4}=0.0025, \quad g_{ZS_3S_4}=0.37. \end{array} \end{aligned}$$
(41)
Fig. 4
figure 4

The \( e^{+} e^{-} \rightarrow Z^{*} \rightarrow S_{i} S_{j}\) cross section for BP C with masses in GeV

Full size image

Points from the C type benchmark scenario, as shown in Fig. 4, have the specific characteristic of seeing all thresholds (nearly) overlapping at low \(\sqrt{s}\) values. The size of the various cross sections is again dictated by the \(ZS_iS_j\) couplings and thus is largest for \(S_1S_2\) and \(S_3S_4\) over any of \(S_1S_3\), \(S_1S_4\), \(S_2S_3\) and \(S_2S_4\), as previously seen already.

4.2 Significance of the signal over the SM background

Linear colliders are perfect machines to test the inert sector of multi-scalar models, e.g. in leptonic or semi-leptonic channels with missing transverse energy. Full analysis is beyond the scope of this work. However, as shown in [49] in the context of the IDM, with the use of multi-variate analysis and proper cut selection it is possible to claim \(5\sigma \) discovery for \(\sqrt{s} = 500\hbox { GeV}\) and integrated luminosity of \(1\hbox { ab}^{-1}\) if the sum of masses of the neutral scalar pair is below 330 GeV [50], which all of our proposed benchmarks fulfil.

With all our scalars lighter than 160 GeV, the final state of the \(e^+e^- \rightarrow Z^* \rightarrow S_i S_j\) process would be missing transverse energy and \(f \bar{f}\),

$$\begin{aligned}&e^+e^- \rightarrow Z^* \rightarrow S_1 S_j \rightarrow S_1 S_1 Z^* \rightarrow S_1 S_1 f \bar{f}, \nonumber \\&e^+e^- \rightarrow Z^* \rightarrow S_i S_j \rightarrow S_1 Z^* S_1 Z^* \rightarrow S_1 S_1 f \bar{f} f \bar{f}, \nonumber \\&\qquad \qquad (i,j=2,3,4). \end{aligned}$$
(42)

The main SM background to such processes is through the channels

(43)

Figure 5 shows the SM cross section for these processes in the energy range relevant to our analysis. Note that this background decreases with increasing energy, whereas the signal is asymptotically flat, as shown in Figs. 2, 3 and 4.

Fig. 5
figure 5

The cross section for the SM background processes to our signal

Full size image

Moreover, a simple calculation yields the cross-section for the aforementioned processes to be

$$\begin{aligned}&\sigma ^{SM}(e^+ e^- \rightarrow W^+W^-) \times \text{ Br }(W^+ \rightarrow l^+ \nu ) \nonumber \\&\qquad \times \, \text{ Br }(W^- \rightarrow l^- \bar{\nu }) \nonumber \\&\quad = 0.09 \; \sigma ^{SM}(e^+ e^- \rightarrow W^+W^-) \lessapprox 1.8 \; \text{[pb] }, \nonumber \\&\sigma ^{SM}(e^+ e^- \rightarrow Z Z) \times \text{ Br }(Z \rightarrow f \bar{f}) \times \text{ Br }(Z \rightarrow \nu \bar{\nu }) \times 2 \nonumber \\&\quad = 0.28\; \sigma ^{SM}(e^+ e^- \rightarrow Z Z) \lessapprox 0.56 \; \text{[pb] }, \nonumber \\&\sigma ^{SM}(e^+ e^- \rightarrow Z h) \times \left[ \text{ Br }(Z \rightarrow f \bar{f}) \times \text{ Br }(h \rightarrow inv.) \right. \nonumber \\&\qquad \left. +\, \text{ Br }(Z \rightarrow \nu \bar{\nu }) \times \text{ Br }(h \rightarrow f \bar{f}) \right] \nonumber \\&\quad = 0.28 \; \sigma ^{SM}(e^+ e^- \rightarrow Z h) \lessapprox 0.112 \; \text{[pb] }, \end{aligned}$$
(44)

which should lead to a non-negligible significance of the signal over the background.

5 Conclusion and outlook

We have studied distinctive signatures of the CP-violating I(2+1)HDM at a future \(e^+e^-\) collider. The off-shell Z boson in the process \(e^+e^-\rightarrow Z^* \rightarrow S_iS_j\) leads to six possible final states involving pairs of dark (or inert) neutral states, \(S_1S_{2,3,4}\), \(S_2S_{3,4}\), \(S_3S_4\), in the CP-violating case, in comparison to four possible final states, \(H_1 A_{1,2}, H_2A_{1,2}\), in the CP-conserving case. We then have provided several BPs, for which we have shown production rates as large as a few picobarns at \(\sqrt{s}\) energies accessible by future electron-positron colliders. The relative distance (in \(\sqrt{s}\)) of the production thresholds of these final states as well as their heights would serve the purpose of separating typical dark scalar mass patterns, of which we benchmarked here three types, each corresponding to different DM dynamics compatible with relic density as well as both direct and indirect searches.

Given the foreseen timescale for the construction and exploitation stage of future \(e^+ e^-\) colliders, we will therefore be able to probe the described I(2+1)HDM benchmark scenarios on time scales of ten to twenty years from now. By that time, we expect an increased sensitivity of DM (in)direct detection experiments and more stringent constraints from a high energy and/or luminosity LHC, so that, by combining information from all these sources, one may be in a position to eventually use extremely collimated and energetically precise electron-positron beams in order to perform a threshold scan able to accurately extract the six rest masses, \(m_{S_i}+m_{S_j}\), in turn leading to a fit to the individual ones, \(m_{S_i}\).

In fact, such a scope offered by future \(e^+e^-\) colliders is complementary to the one that will be afforded by, e.g., a XENONnT upgrade. Furthermore, the former experiments are very useful for testing the I(2+1)HDM, as their measurements do not depend on the small Higgs-DM coupling, unlike the latter. Finally, it is worth noting that XENONnT (and other direct detection experiments) are sensitive only to the mass and couplings of the DM particle and will provide no information about other unstable particles from the dark sector. Therefore, future \(e^+e^-\) colliders will be essential to probe other inert particle masses than the DM candidate one, all of which will hardly be accessible at the LHC [46], for the simple reason that, at an ILC [47] or FCC-ee [48], it is the highly controlled initial state that enable access to these while, at the LHC, this happens through final states always containing missing energy.