Exploring collider aspects of a neutrinophilic Higgs doublet model in multilepton channels

We consider a neutrinophilic Higgs scenario where the Standard Model is extended by one additional Higgs doublet and three generations of singlet right-handed Majorana neutrinos. Light neutrino masses are generated through mixing with the heavy neutrinos via Type-I seesaw mechanism when the neutrinophilic Higgs gets a vacuum expectation value (VEV). The Dirac neutrino Yukawa coupling in this scenario can be sizable compared to those in the canonical Type-I seesaw mechanism owing to the small neutrinophilic Higgs VEV giving rise to interesting phenomenological consequences. We have explored various signal regions likely to provide a hint of such a scenario at the LHC as well as at future $e^+e^-$ colliders. We have also highlighted the consequences of light neutrino mass hierarchies in collider phenomenology that can complement the findings of neutrino oscillation experiments.


I. INTRODUCTION
The discovery of the 125 GeV Higgs boson [1,2] has been a remarkable achievement of the Large Hadron Collider (LHC). This has provided us a closure regarding the predictions of the Standard Model (SM). While our quest towards understanding the physics beyond the Standard Model (BSM) continues, the 13 TeV run of the LHC is expected to make a big impact both in terms of higher energy reach and better precision by accumulating huge amount of data at large luminosity. The enigma of non-zero neutrino mass has pushed the theorists as well as experimentalists to develop new theories and experimental techniques in order to establish the right theoretical pathway towards unveiling the true nature of neutrino mass generation. The neutrino oscillation experiments have established the fact that at least two of the three light neutrinos are massive, and that they have sizable mixing among themselves (for a review, see [3]). The SM, lacking any right-handed neutrinos, is unable to account for these phenomena. This has led to a plethora of scenarios leading to neutrino mass generation [4][5][6][7][8][9][10][11][12]. As the resulting neutrino mass eigenstates may be either Dirac or Majorana type, both scenarios have potentially unique signatures [13][14][15][16][17][18][19][20][21][22][23] in the collider experiments. The LHC collaborations have put significant effort to extract any possible information about such scenarios from the accumulated data and the null results so far have only been able to constrain the parameter space of various neutrino mass models [24][25][26][27][28][29].
In the post Higgs discovery LHC era, the true nature of the scalar sector remains another vital area of interest. The natural question that arises is whether the 125 GeV Higgs is the only scalar as predicted by the SM or other exotic scalars exist alongside, as predicted by various BSM theories including some of the neutrino mass models [10,30]. The measurements of couplings of the 125 GeV Higgs with known SM particles have so far been consistent with the SM predictions [31]. Thus, even if this Higgs boson were indeed part of a larger scalar sector, its mixing with the other states would be small. There are still enough uncertainties in these measurements to allow new exotic scalar multiplets. Unless the LHC observes some hint of a new scalar, our only hope lies in the precision measurements of the Higgs couplings in order to constrain the BSM physics scenarios. Meanwhile, there has been a long term interest in the simplest two-Higgs doublet models (2HDM) (for a review, see [32]) which are also strongly motivated by supersymmetric scenarios. A two-Higgs doublet model predicts the presence of two CP-even, one CP-odd and two charged Higgses, one of the CP-even Higgs states being the 125 GeV Higgs boson. Despite the presence of these additional scalar states, the mixing between the two doublets can be arranged so that the other scalars are practically decoupled from the SM Higgs. In such cases, the interaction of the SM-like Higgs with the exotic scalars may be so suppressed that any hint of such interactions can be very hard to pick up even with the precision measurements at the LHC.
The hope of finding these scalars, therefore, lies in their direct search. While the increasing center-of-mass energy at the LHC can probe heavier exotic particles, extracting any new physics information from the tremendous amount of collected data also faces the increasing challenge of tackling the QCD background. Hence looking for lepton-enriched final states is understandably efficient in suppressing the SM background contributions and probing new physics scenarios which can potentially give rise to lepton-rich final states.
In this work, we consider a 2HDM where the additional Higgs doublet has an odd Z 2 symmetry charge opposite to all the SM particles, preventing it from interacting directly with the leptons and quarks. One can additionally incorporate right-handed neutrinos in the model with similar transformation property under Z 2 symmetry as the new Higgs doublet.
One can thus generate Dirac neutrino mass terms when the Z 2 breaks spontaneously and the new Higgs doublet gets a vacuum expectation value (VEV). This class of models, known as neutrinophilic Higgs doublet models (νHDM) have been proposed long ago [33][34][35] and the relevant phenomenology has been studied quite extensively [36][37][38][39][40][41][42][43][44]. In principle, one can also generate Majorana neutrino mass terms in such a scenario, since a Majorana mass term for the additional right-handed neutrinos does not break the Z 2 symmetry but breaks the accidental lepton number symmetry by two units (∆L = 2). Such neutrino mass generation mechanism looks very similar to the Type-I seesaw [4][5][6][7] case, save for the fact that one uses the neutrinophilic Higgs VEV instead of electroweak VEV in order to generate the lightheavy neutrino mixing. The advantage of having the additional Higgs doublet to generate non-zero neutrino masses is that the additional VEV can be very small 1 in order to counter the smallness of the light neutrino masses which would otherwise be fit with a very small Dirac neutrino Yukawa coupling that has no significant collider phenomenological aspects.
Depending on whether the non-zero neutrinos are Dirac or Majorana type, the collider signals of a νHDM scenario can be very different. When Majorana neutrinos exist, smoking gun signal would be lepton number violating final states. In this work, instead of looking for direct heavy neutrino production, we have considered the production of the neutrinophilic charged Higgs (H ± ) and explored its various possible decay modes. There are some earlier studies on the charged Higgs in similar scenarios emphasising its decay into a charged lepton and a heavy neutrino in the process [44]. We show that even cleaner signals can be obtained using this decay mode with higher lepton multiplicity where the SM background is practically non-existent. We also show that sizable signal event rates can be obtained with other possible decay modes of the H ± , which can serve as complementary channels in probing a νHDM-like scenario. We perform our analysis using the 13 TeV LHC as well as an e + e − collider with 1 TeV center-of-mass energy. In the process, one can extract information on the neutrino sector parameters also. We show that a very clean indication of the neutrino mass hierarchy can be obtained from the multiplicity of the charged leptons in the final state even after a rigorous collider simulation. Such information can be very useful in complementing the neutrino oscillation experiments.

II. MODEL
In the νHDM model, the particle content of the SM is extended by one additional Higgs doublet (φ ν ) and three generations of SM gauge singlet right-handed neutrinos (N ). A discrete Z 2 symmetry is introduced, under which both φ ν and N i , i = 1, 2, 3, are odd while all the SM fields are even. The most general scalar potential involving the two Higgs doublets is given by where a non-zero m 3 explicitly breaks the Z 2 symmetry in the model. In the absence of this term, the Z 2 symmetry can be broken spontaneously by a vacuum expectation value (VEV) v ν of the field φ ν , while the standard electroweak symmetry is broken when φ acquires a VEV, v φ .
Let us first discuss a framework, where m 3 = 0, i.e. Z 2 symmetry is broken only spontaneously in order to generate light neutrino masses and mixing. The model is constrained by sterile neutrino searches, effective number of neutrinos and amount of 4 He required in Due to an instability of right-handed neutrinos N i induced by their mixing to left-handed neutrinos, the mixing strength between ν , = e, µ, τ , and N i , that is, |U i |, can be probed by sterile neutrino searches. In semileptonic meson decays, N i are produced, and can subsequently decay to charged leptons and mesons. Present constraints on |U ei | 2 and |U µi | 2 allow a region where their magnitude is of order 10 −10 to 10 −6 , assuming M N i < 2 GeV [29]. For tau-sterile mixing, |U τ i | 2 10 −4 , assuming M N i < 0.3 GeV.
In νHDM, however, we found the model favoring even lower values of active-sterile mixing, of order |U i | 2 ∼ 10 −18 to 10 −12 , at M N i = 1 GeV, and even lower for higher Majorana neutrino masses (see Fig. 1). The largest and smallest active-sterile mixings are driven by U τ 3 and U τ 1 elements. Therefore all the active-sterile mixing elements fall between them: |U τ 1 | < |U i | < |U τ 3 |. The matrix elements are proportional to M −1/2 N , therefore M N with some constants C 1 and C 2 . They are deduced from Fig. 1, having values C 1 = 4.34 · 10 −9 MeV 1/2 and C 2 = 2.34 · 10 −6 MeV 1/2 . The matrix elements then belong to the following interval: In addition, the constraints for |U i | were derived from assumption that the branching ratios for N i decay are dominant. This is not applicable for νHDM, since then the decay modes of right-handed neutrinos are dominated by decays to invisible particles.
As the model is unconstrained by semileptonic and leptonic decay modes, the lower bound for M N i arises from BBN. In the early universe the right-handed neutrinos must be heavy enough to fall off from the thermal equilibrium before BBN. This is due to the latest results for effective number of neutrinos (N ν = 3.15 ± 0.23) by PLANCK [46], which forbids large interference from right-handed neutrinos. This leads to a constraint M N i 100 MeV.
In addition neutrinophilic VEV v ν is constrained from both above and below. Ultralight VEV is forbidden by astrophysical constraints: v ν O(eV) [47,48]. On the other hand, the surface energy density associated with the domain wall arising from discrete Z 2 symmetry breaking is η ∼ v 3 ν [49]. The effect of these domain walls to the temperature anisotropies of CMB is where G is Newton's gravitational constant, H 0 is Hubble constant and we have assumed [50]. Since the observed temperature anisotropies by PLANCK are ∼ 10 −5 , the birth of a domain wall will not contradict cosmological data if the VEV is small. If we require the contribution to CMB temperature anisotropies not to exceed the experimental limit, together with the astrophysical constraints, we get In order to apply perturbative theory to νHDM, the absolute values of the elements of the light neutrino Yukawa coupling matrices must be √ 4π at most. We performed a global fit to available neutrino oscillation data to calculate the matrix elements, assuming normal Blue-shaded region denoted 'CMB' is excluded due to restrictions of CMB temperature anisotropies induced by domain walls. The available parameter space is restricted also from BBN requirement neutrino mass ordering, higher θ 23 octant and no CP violation. We found the dependence of the largest Yukawa coupling of v ν and M N to be The dependence is illustrated in Fig. 2.
The breaking of Z 2 symmetry is necessary in order to generate light neutrino masses within the framework of this model by means of their mixing with heavy right-handed neutrinos. One can add the following Yukawa interaction and Majorana neutrino mass terms to the Lagrangian while keeping the Z 2 parity unbroken: where, m ij R represents the Majorana mass terms corresponding to the right-handed neutrinos. Once φ ν acquires a VEV, the Yukawa term gives rise to Dirac neutrino mass terms, m ij D = v ν × y ij ν . The physical Higgs sector now consists of two neutral CP-even (h, H ν ), one neutral CP-odd (A ν ) and the charged Higgs (H ± ν ) 2 . In the case when m 3 = 0, the physical mass eigenvalues at tree level are given by: . v ν being small, terms proportional to v n ν (where n > 2) have been neglected. Note that the mixing angle between the SM and neutrinophilic Higgs states are proportional to the ratio vν v φ and can be safely neglected since we assume v ν v φ . Under this circumstance, the CP-even neutrinophilic Higgs (H ν ) is always light and the heavy neutrino almost always decay into H ν and a light neutrino resulting in an opposite-sign dilepton signal for a charged Higgs pair production channel [40]. However, if the explicit symmetry breaking term is present in the Lagrangian, i,e, m 3 = 0, the mass eigenvalues are given by: Now the neutrinophilic CP-even Higgs can be heavy depending on our choice of m 3 , m Hν A heavy H ν and (or) A ν opens up the possibility of a cascade decay via heavy neutrinos resulting in multi-lepton signals of such a scenario that we intend to explore. In the limit m 3 → 0, the symmetry of the theory is enhanced. Thus, m 3 can be assumed to be naturally small. Besides, a large m 3 can also give rise to significant mixing between the two Higgs doublets, which is strictly constrained from the present Higgs data.

A. Neutrino Mass Generation
The neutrino oscillation data [3,51,52] indicates that at least two of the three light neutrinos have non-zero mass. One of the most natural ways to generate tiny neutrino mass 2 We have assumed the scalar potential to be CP invariant.
is via seesaw mechanism [4][5][6][7][8][9][10][11][12]. In νHDM the mechanism is very similar to that of Type-I seesaw [4][5][6][7]. The mixing between light and heavy neutrinos is introduced via the term y νL φ ν N in the aftermath of symmetry breaking, when φ ν gets a VEV. In the basis {ν, N } the 6 × 6 neutrino mass matrix looks like where m D = y ν v ν . The light effective 3 × 3 neutrino mass matrix in the approximation m D << m R is given by The above equation looks exactly similar to what we obtain in canonical Type-I seesaw scenario. The only difference is that in the present framework v ν can be quite small and as a result one can have larger y ν compared to the canonical Type-I seesaw scenario, thus making this model phenomenologically more interesting. In order to fit the oscillation data, one also needs to account for the mixing among the three light neutrino states constrained by the PMNS matrix. One can rewrite M ν in Eq. (10) as where m diag ν is the diagonal light 3 × 3 neutrino mass matrix and U is the PMNS mixing matrix. In order to produce proper mixing satisfying the experimental bounds on the PMNS matrix elements, one of the matrices, m D or m R , has to be off-diagonal. Here we choose to keep m R diagonal and fit the PMNS matrix via an off-diagonal m D . Thus y ν is obtained using Casas-Ibarra parameterization [53] where R can be any orthogonal matrix and complex provided R T R = 1. For simplicity, we have chosen R to be an identity matrix.
Thus with correct choices of the parameters m D and m R , Eq. (9) is capable of explaining the neutrino oscillation data at the tree level itself. There is a potential source of large correction [54,55] [37,[56][57][58]. As can be seen both from Eq. (7) and (8), the mass splitting between these two states is driven by the parameter λ 5 which is therefore set equal to zero throughout this work.

III. CONSTRAINTS AND BENCHMARK POINTS
Constraints on the charged Higgs mass and its couplings may arise from direct collider search results, neutrino oscillation data and lepton flavor violating decay branching ratios.
The LHC collaborations have looked for signatures of exotic scalars in various channels and put bounds on the charged Higgs mass in the range 300 -1000 GeV provided it can decay only into a top and a bottom quark [59][60][61][62]. However, in our present scenario, the charged Higgs, being a neutrinophilic one, does not couple to the quarks. In such scenarios, there are no direct search constraints on m H ± ν . In principle, the constraints derived from slepton searches at the LHC can be reinterpreted to put bounds on the neutrinophilic charged Higgs masses although only in the massless limit of the lightest neutralino. Two body decay of the sleptons into a charged lepton and lightest neutralino gives rise to a dilepton signal which can be relevant for the the present scenario. Existing data excludes slepton masses upto 450 GeV in presence of massless neutralino [63,64]. However, one always obtains same-flavor-opposite-sign(SFOS) lepton pairs from such slepton pair production processes.
The signal requirement also demands a jet veto in the central region alongside the SFOS lepton pair for such analyses. In the present scenario, largest event rate in such a signal region can be obtained when H ± ν decays into a charged lepton and a heavy neutrino. Heavy neutrino further decays into a light neutrino and Z-boson which further decays invisibly.
Clearly, the resulting signal cross-section is rendered small due to branching suppressions.
Demand of SFOS lepton pairs makes this cross-section even smaller 3 . Thus, the existing slepton mass limit when reinterpreted for m H ± ν proves to be much weaker. Its couplings with the heavy neutrinos on the other hand, can be constrained from neutrino oscillation data and lepton flavor violating decay branching ratios [44]. As mentioned in section II A, we have used off-diagonal m D while fitting the PMNS matrix. These off-diagonal entries are severely constrained from LFV decay branching ratio constraints [65][66][67][68][69][70][71]. These constraints are also reflected upon our choice of the neutrinophilic Higgs VEV, v ν . It has been observed and also verified by us that v ν can be ∼ 10 −2 GeV [44] at the smallest, if the neutrino oscillation data and the LFV constraints are to be satisfied simultaneously, the most stringent constraint arising from the non-observation of BR(µ → eγ) [65,66]. This constraint puts the spontaneously breaking Z 2 scenario in jeopardy. As evident from Fig. 2, such a choice of v ν is clearly ruled out from restrictions on CMB temperature anisotropies induced by domain walls. However, if the Z 2 symmetry is broken explicitly, this domain wall problem can be averted. Hence for this work, we choose to work with the m 3 = 0 scenario only.

A. Charged Higgs branching ratios and pair production cross-section
The possible decay modes of the neutrinophilic charged Higgs (H ± ν ) in our present scenario The relevant interaction vertices are given in the Appendix. Depending on the mass hierarchy of H ± ν , H ν (A ν ) and N , and the choice of neutrino mass hierarchy one (or two) of these decay modes determine the event rates of the different possible final states at the collider. Note that the branching ratios of the decays into the neutral CP-even and CP-odd Higgs states are always the same since they are mass degenerate by our choice of the parameters. These two decay modes dominate over the heavy neutrino decay modes always, if the mass difference, ∆m = m H ± ν −m Hν (m Aν ) is larger than that of the W-boson mass, m W . This is an artefact of the small Dirac neutrino Yukawa parameters, which are otherwise constrained by neutrino oscillation data and the non-observation of LFV decays. The y ν being smaller by orders of magnitude from the competitive gauge coupling, a large branching ratio into the N ± or N τ ± decay modes are not ensured even if ∆m < m W . In spite of the additional phase space suppression, three-body decays of H ± ν via off-shell W -decay, dominate over these two-body modes unless ∆m m W .
This behavior is depicted in Fig. 3 where the competitive nature of BR(H ± ν → N ± ) and BR(H ± ν → H ν (A ν ) ν), where = e, µ, is clearly visible through the distributions of the starred and circular points respectively. BR(H ± ν → N ± ) overtakes the three-body decay branching ratio only if ∆m < 50 GeV. For ∆m > m W , BR(H ± ν → H ν (A ν )W ± (→ ± ν)) takes over and remains the only dominant decay mode.
In Fig. 4, we have shown variation of the H ± ν production cross-section at the LHC and an e + e − collider. The figure on the left shows the variation of the cross-sections as a function of m H ± ν at the 13 TeV LHC and different center-of-mass energies (500 GeV, 1 TeV and 3 TeV) at an e + e − collider. Note that, at the LHC, the H ± ν production channels include pp → H ± ν H ∓ ν , pp → H ± ν H ν and pp → H ± ν A ν while for the e + e − collider, pair production is the only viable option. Since we have assumed m Hν = m Aν for our study, the crosssections of the above mentioned second and third production channels exactly equal. Hence we have shown their combined cross-section in the figure and evidently, it dominates over the pair production cross-section throughout the entire charged Higgs mass range. However, both these cross-sections fall rapidly with increasing mass. On the other hand, at an e + e − collider the cross-section falls far less rapidly implying the fact that such a collider will be collider is likely to be much cleaner in terms of the SM background contributions. In this work, we have taken into account all the aforementioned production channels for LHC and just the pair production for the e + e − collider analysis.

B. Choice of benchmark points
We now proceed to choose some benchmark points representing the different interesting features of the present scenario for further collider studies. As discussed earlier, one can obtain different possible final states depending upon the mass hierarchies of H ± ν , H ν (A ν ) and N . Since we also aim to correlate the light neutrino mass hierarchy with the multiplicity of different lepton flavor final states, we will study cases in which at least one of the heavy neutrinos is lighter than the neutrinophilic Higgs states so that it can appear in the cascade.
In Table I below we present the input parameters, relevant masses and the resulting y ν for the four benchmark points of our choice. We have incorporated the complete model in SARAH [72][73][74][75][76], and subsequently imported in SPheno [77,78] in order to perform the analytical and numerical computation of the masses and mixings of the particles, their branching ratios and other relevant constraints. See Appendix for LFV constraints for our benchmarks.  The four benchmark points are chosen such that all the dominant decay modes of the neutrinophilic Higgs and the heavy neutrinos are highlighted by different mass hierarchies.
The relevant branching ratios are shown in Table II. The two most dominant decay modes of H ± ν are N , where = e, µ, τ , and H ν (A ν )W ± . The first decay mode is driven by the Dirac neutrino Yukawa couplings, y ν , whereas the second one is driven by gauge couplings. As discussed above, the elements of y ν are already constrained from the neutrino oscillation data as well as from the LFV constraints, and thus are in general weaker than the competitive gauge coupling. Hence, if the mass splittings among the neutral and charged neutrinophilic Higgs and the heavy neutrino states are such that both N and H ν (A ν )W ± decay modes are kinematically accessible for H ± ν , the gauge boson associated one becomes its only relevant decay mode. However, if at least one of the heavy neutrinos is lighter than the H ± ν and the H ν (A ν ) states are almost degenerate to it, then the decay via heavy neutrinos becomes important. The latter scenario is highlighted in BP1 and BP3 while BP2 represents the former scenario. BP4, on the other hand, highlights the situation where the two-body mode N competes with the three-body decay into H ν (or A ν ) alongside an off-shell W -boson.
However, ∆m being on the larger side, the three-body decay dominates as discussed earlier in Section III A. The heavy neutrinos (N ) in this scenario can decay either via the SM gauge bosons (W ± , Z) or the different Higgs states. Note that decays of N into W ± , Z and h can only occour through their mixing with the light neutrinos which are suppressed in the present scenario. Hence, these decay modes become relevant for N only if the neutrinophilic Higgs states are kinematically inaccessible to it. The choice of neutrino mass hierarchy clearly reflects in the branching ratios of both H ± ν and N and is also expected to be reflected in the final event rates of the multi-lepton signals we intend to explore.

A. Identifying signal regions
In the context of LHC, we aim to study cleaner multi-lepton channels with no tagged jets in the final state. The possible final states that we probe in the present context are , and same-sign trilepton (SS3 ) + E / T + X (SR3), where X represents everything else (jets, photons or leptons) 4 in the final states.
As mentioned earlier, the various branching ratios of H ± ν and hence the final signal event rates depend on the mass difference factor ∆m. Thus, it is interesting to study how the signal rates vary depending on ∆m which in turn can also provide an indirect hint about the masses of H ν and (or) A ν .
In Fig. 6 we have shown the variation of the cross-sections corresponding to the three signal regions mentioned above as a function of ∆m with color-coded m H ± ν . Note that these cross-sections are theoretical estimates obtained after combining contributions from all the three relevant production modes of H ± ν at the LHC prior to detector simulation and do not include the cut efficiencies. The two rows of the figures correspond to normal and inverted hierarchy of neutrino masses respectively. While SR1 only receives contribution from pairproduction, both SR2 and SR3 are enriched with contributions from pair production as well as associated production of the H ± ν . Most of the signal events corresponding to SR1 and SR2 are expected to arise from H ± ν decay into a charged lepton and a heavy neutrino followed by the heavy neutrino decay into a charged lepton and W . Depending on the leptonic or hadronic decays of the W -bosons, one can obtain various lepton multiplicities as represented by these signal regions. The signal cross-sections are largest when H ± ν and H ν (A ν ) are mass degenerate for any given m H ± ν and they drop with increasing m H ± ν and ∆m. SR3, on the other hand, receives more contribution when H ± ν decays into H ν (or A ν ) along with an on-shell or off-shell W -boson. Moreover, in the case of pair production, same-sign leptons can not be obtained if both the H ± ν decays via ± N resulting the crosssection of SR3 being smaller with smaller ∆m. However, the associated production channels contribute dominantly to this signal region throughout the whole range of ∆m. As evident from Fig. 6, SR3 is the most favorable channel to look for such scenarios. In general, the inverted hierarchy of the light neutrino masses is expected to generate more multi-leptonic events owing to the larger branching ratio, BR(H ± ν → eN ) as reflected by the plots on the bottom row.
In a similar way, we now proceed to choose some signal regions for our analyses in the context of an e + e − collider. The possible final states we probe in this context are ≥ 5 +E /+X (SR4), ≥ 4 + ≥ 2 − jet + E / + X (SR5), and SS3 + E / + X (SR6) 5 . The corresponding signal rates are showcased as a function of ∆m in Fig. 7. Trends of the distributions are similar to what we obtained for the LHC case. However, the difference in the production cross-section is manifested by the signal cross-sections indicating a larger event rate at the LHC for similar final states at low m H ± ν region. The rapid fall in production cross-section with increasing m H ± ν at the LHC makes it less relevant for heavier charged Higgs masses.
An e + e − collider can be more effective provided the center-of-mass energy is large enough for the production. Here, the signal rates drop alarmingly close to m H ± ν ∼ 500 GeV due to the choice of center-of-mass energy as 1 TeV.

B. Analysis
In order to carry out the simulation, events were generated at the parton level using MadGraph5 [80,81] with nn23lo1 parton distribution function [82,83] and the default dynamic factorisation and renormalisation scales [84]. We have used PYTHIA [85] for the subsequent decay of the particles, showering and hadronization. After that the events are passed through Delphes [86][87][88] for detector simulation. Jets have been reconstructed using anti-kT algorithm via FastJet [89,90]. The b-jet and τ -jet tagging efficiencies as well as the mistagging efficiencies of the light jets as b or τ -jet have been incorporated according to the latest ATLAS studies in this regard [91].

Primary selection criteria
We have applied the following cuts (C0) on the jets, leptons and photons in order to identify them as final state particles: • All the charged leptons are selected with a transverse momentum threshold p T > 10 GeV and in the pseudo-rapidity window |η| < 2.5.
• All the jets including b-jets and τ -jets must have p j T > 20 GeV and |η| j < 2.5.
• We demand ∆R ij > 0.4 between all possible pairs of the final state particles to make sure they are well separated.
As discussed in Section III, the choice of neutrino mass hierarchy affects the branching ratios of the neutrinophilic Higgs as well as the heavy neutrinos in certain flavor specific decay modes. Thus, the hierarchical effect is reflected by the abundance of certain flavor of leptons in the signal events. As we have seen, one would expect less abundance of electrons in the final states for normal hierarchy scenario compared to that for inverted hierarchy.
This feature is evident in Fig. 8 which shows the electron multiplicity in the final state with at least four leptons for BP1 in normal as well as inverted hierarchy scenarios. Such lepton multiplicity distributions can thus provide indirect probe of existing neutrino mass hierarchy.

C. Results@LHC13
In the context of LHC, we have studied the final states corresponding to SR1, SR2 and SR3 as defined in Section IV A. Although the choice of our signal regions ensure small or no SM background, we have checked the relevant production channels, tt, ttZ(γ * ), ttW , W W Z, W ZZ, ZZZ and ZZ + jets nevertheless in this regard. In Table III   As expected SR1 has the smallest event rate owing to its large lepton multiplicity, but with negligible SM background. Thus it can be a very clean signal but only if the charged Higgs mass is on the lighter side, as in BP1, and at least one of the heavy neutrinos is lighter than the charged Higgs. The situation, however, worsens considerably with increasing charged Higgs mass, as indicated by BP3. SR2 has a much better event rate and can probe BP1 and BP3 at much lower luminosity than SR1. As the numbers in Table III indicate, the inverted hierarchy scenario for BP1 can be probed with a 3σ statistical significance at an integrated luminosity of ∼ 30 fb −1 in both SR2 and SR3, i.e, if these signal regions are studied, m H ± ν in this mass range can be probed and possibly be excluded with the LHC data already accumulated. For the corresponding normal hierarchy case however, for the same benchmark point, one needs L ∼ 50 fb −1 for similar discovery significance in SR3. BP3 requires an integrated luminosity of ∼ 100 fb −1 (for inverted hierarchy) or more. For the benchmark points like BP4 and BP2, the decay H ± ν → N is either suppressed or absent altogether. Thus for such points SR1 ceases to be a viable signal region while SR2 is relevant only at large luminosities. In this case SR3 turns out to be the most viable signal region. In this signal region, to achieve 3σ statistical significance in the inverted hierarchy case of BP4 and BP2 one requires L ∼ 100 fb −1 and ∼ 200 fb −1 respectively. For all the benchmark points, the choice of light neutrino mass hierarchy is clearly manifested through the different signal event rates. Evidently, with multi-leptonic final states, an inverted hierarchy scenario is more likely to be probed at lower luminosities at the LHC.
In Table IV below Table IV correspond to L = 100 fb −1 , the inverted hierarchy scenarios in BP1 and BP3 can be probed with a statistical significance of 3σ at a much lower luminosity (∼ 10 fb −1 ). Note the improved event rates in signal regions SR5 and SR6 despite of the smaller H ± ν pair production cross-section in BP3 over those of BP1. This is a consequence of increased hadronic branching ratio of N and improved cut efficiency due to the larger mass gap between H ± ν and N in BP3. Even BP2 which can be probed at the LHC only at very high luminosity can be probed here at around L = 50 fb −1 with similar statistical significance via SR5 which turns out to be the most favored signal in general for all the benchmark points. The overall signal rate is relatively weaker in SR6 due to better lepton tagging efficiency at a lepton collider, which results in smaller number of events with exactly three same-sign leptons as demanded.  We have also highlighted the interesting roleplay of the light neutrino mass hierarchy.
Whether the neutrinos follow normal or inverted hierarchy, is likely to be manifested via multiplicity of different flavored leptons in the final state. Thus such a finding at the collider experiments can complement the neutrino oscillation experiments which are yet to ascertain the correct mass hierarchy of the three light neutrinos.
The fact that the charged Higgs pair production cross-section falls quite rapidly at the LHC with increasing mass, led us to perform a comparative study between the LHC and a future e + e − machine in order to probe such scenarios. We observed that although LHC is quite efficient to probe light charged Higgs masses, an e + e − collider will be able to probe a much larger parameter space with heavier states.
The work of SKR was partially supported by funding available from the Department of Atomic Energy, Government of India, for the Regional Centre for Accelerator-based Particle Physics (RECAPP), Harish-Chandra Research Institute.

VII. APPENDIX
A. Relevant interaction vertices where U ij L , U ij N and Z ij are the charged lepton, neutrino and CP-even neutral Higgs mixing matrices, the bases of the mass matrices being {e, µ, τ }, {ν e , ν µ , ν τ , N c e , N c µ , N c τ } and {φ, φ ν } respectively. Note that, U ij L is a diagonal matrix. In Table V we have shown the obtained branching ratios for various lepton flavor violating processes corresponding to the four benchmark points. BR(µ → eγ) is projected to be probed experimentally up to 6.0 × 10 −14 in near future [66]. As indicated by the numbers, the obtained branching ratios for this process are at least one order of magnitude smaller for our benchmark points. The rest of these obtained LFV branching ratios are several orders of magnitude below the present experimental sensitivity in the respective channels.