Probing Electroweakly Interacting Massive Particles with Drell-Yan Process at 100 TeV Hadron Colliders

There are many models beyond the standard model which include electroweakly interacting massive particles (EWIMPs), often in the context of the dark matter. In this paper, we study the indirect search of EWIMPs using a precise measurement of the Drell-Yan cross sections at future $100\,{\rm TeV}$ hadron colliders. It is revealed that this search strategy is suitable in particular for Higgsino and that the Higgsino mass up to about $1.3\,{\rm TeV}$ will be covered at $95\,\%$ C.L. irrespective of the chargino and neutralino mass difference. We also show that the study of the Drell-Yan process provides important and independent information about every kind of EWIMP in addition to Higgsino.


Introduction
There are many models that extend the standard model (SM) by introducing electroweakly interacting massive particles (EWIMPs). One motivation to introduce them is the existence of the dark matter (DM) in our universe. For example, an electroweakino in the minimal supersymmetric (SUSY) extension of the SM can be the lightest supersymmetric particles (LSP) and are natural DM candidate. In particular, there are well motivated scenarios where the LSP is Higgsino or wino, which transform as doublet and triplet under the weak SU (2) L gauge symmetry, respectively; light Higgsino is preferred to reduce the amount of the finetuning of the electroweak scale as in the "natural SUSY" set up [1][2][3][4], while the so-called "mini split" spectrum [5][6][7][8][9][10] with anomaly mediation [11,12] makes wino LSP. Another example is the "minimal dark matter" (MDM) scenario [13][14][15], in which an additional EWIMP is introduced, whose stability is ensured due to the large SU (2) L charge assignment. In particular, 5-plet fermion and 7-plet scalar with hypercharge zero are good DM candidates that escape from the DM search experiments so far.
In order to search for EWIMPs, several different approaches are adopted. One approach is the direct production at the large hadron collider (LHC). The main strategy is to use the disappearing charged track, which indicates a long-lived charged particle, the charged components of the EWIMP in the current case. Both ATLAS and CMS collaborations announced a result of this method with the data of √ s = 13 TeV LHC [16][17][18]. The current lower bound on the mass of the pure Higgsino-like (wino-like) state is 152 (460) GeV at 95 % C.L. We can obtain a similar bound for the MDM using the same method [19]. In this method, however, the bound strongly depends on the lifetime of the charged component, which is sensitive to the mass difference between the charged and the neutral components.
In particular, it is often the case in the SUSY model that the Higgsino-like LSP and its charged counterpart possess a non-negligible fraction of wino, which significantly enhances the mass difference compared to the pure Higgsino case. In such a case, the lifetime of the charged component is extremely short, making the disappearing track method challenging. There is another option called mono-X search to search for a new physics signal in general. However, the corresponding bound is usually very weak due to the large SM background and no bound is imposed on Higgsino at √ s = 14 TeV LHC [20]. Another way of probing EWIMPs is DM search experiments, assuming that the EWIMPs are the dominant component of the DM. Firstly, there exist several direct detection experiments that utilize a scattering between the DM and the nucleon [21][22][23]. Wino is the most promising target of these experiments, whose spin independent scattering cross section with nucleon is almost mass independently given as σ SI 2.3 × 10 −47 cm 2 [24][25][26][27][28], which is still an order of magnitude below the current experimental limit. It is unlikely to detect Higgsino DM in this way since its small SU (2) L charge makes scattering cross section comparable to or below the neutrino floor [25]. The detection of MDM may also be difficult [29] since its possibly larger mass of O(10) TeV weakens the sensitivity of direct detection experiments.
Secondly, a lot of effort is devoted to detecting cosmic rays resulting from DM annihilation, namely the DM indirect detection [30][31][32][33]. Although the results suffer from some astrophysical uncertainties, they have already excluded, for example, wino with mass less than 400 GeV and also around 2 TeV [34]. On the other hand, the corresponding Higgsino bound is weaker and it has been probed only up to 350 GeV [35] again due to the smallness of its SU (2) L charge. For the MDM, 5-plet fermion is analyzed as an example in [36] and the mass less than 2 TeV and several narrow regions are excluded. Note again that the EWIMPs must be the dominant component of the DM for these approaches to be sensitive.
In order to acquire somewhat stronger bounds on the EWIMPs, independent of the decay product or the lifetime and whether they are the dominant component of the DM or not, it has been discussed that indirect search of EWIMPs at collider experiments is useful [37][38][39][40][41]. It utilizes the EWIMP loop effect on various observables. We pursue this possibility and study the effects of EWIMPs on the oblique correction to the electroweak gauge bosons. In particular, in this paper, we study the prospect of the indirect search method at future 100 TeV hadron colliders such as FCC-hh at CERN [42][43][44] or SppC in China [45,46]. We concentrate on the Drell-Yan process that has two leptons in the final state since it provides a very clean signal without any hadronic jets at least from the final state particles. We will show that it provides a comparable or better experimental reach for Higgsino compared to the direct production search of EWIMPs at future colliders [47][48][49]. This method also provides independent and additional information about wino and the MDM.
The rest of the paper is organized as follows. In Sec. 2, we summarize the systematic way of calculating EWIMP effect on the Drell-Yan process. In Sec. 3, we describe our statistical procedure to determine the bounds and show the corresponding results. We also briefly comment on the mass determination of EWIMPs. In Sec. 4, we summarize our results.

EWIMP effect on the Drell-Yan process
We assume that all the beyond SM particles except for EWIMPs are much heavier than O(1) TeV. Then, new physics affects the Drell-Yan process only through the vacuum polarization of the electroweak gauge bosons with EWIMPs in the loop. After integrating out EWIMPs, #1 this effect can be expressed as where L SM is the SM Lagrangian, D is a covariant derivative acting on W aµν , m is the EWIMP mass, #2 and W a µν and B µν are field strength associated with the SU (2) L and U (1) Y #1 Contrary to the usual effective field theory approach, our prescription can also be applied to the energy scale Q m since we do not perform a series expansion of the loop function f in Eq. (1). #2 Here we neglect a small mass splitting among SU (2) L n-plet.
gauge group, respectively. The function f (x) is a loop function defined as where the first (second) line corresponds to a fermionic (scalar) EWIMP, respectively. The coefficients C 1 and C 2 for an n-plet EWIMP with hypercharge Y are given by where g and g are the U (1) Y and SU (2) L gauge couplings, and κ = 1/2, 1, 4, 8 for a real scalar, a complex scalar, a Weyl or a Majorana fermion, and a Dirac fermion, respectively. In particular, (C 1 , C 2 ) = (g 2 /2, g 2 /2) for Higgsino and (C 1 , C 2 ) = (0, g 2 ) for wino. The parton level matrix element for the Drell-Yan process q(p)q(p ) → − (k) + (k ) at the leading order (LO) within the SM is where m and m V are the final state lepton invariant mass and the electroweak gauge boson mass, respectively. Here, {Γ } with e = gs W and g Z = g/c W , where s W ≡ sin θ W , c W ≡ cos θ W with θ W being the weak mixing angle. The parameters Q f , v f , a f for the up-type quarks, the down-type quarks, and the leptons are given by (1), we can also calculate the leading contribution to the matrix element from EWIMPs as where In order to calculate the differential cross section dσ/dm , we need to integrate over the initial parton energy distributions and the final state phase space. Let dL ab /dm be a luminosity function for a fixed m , where a, b denote species of initial partons, √ s is a center of mass energy of the proton collision ( √ s = 100 TeV in our case), and f a (x) is a parton distribution function (PDF) of the given parton a. We also define dΠ LIPS as a Lorentz invariant phase space for the two body final state. Then, with Eqs. (5) and (6), we obtain the differential cross section as Here, dσ SM /dm is the SM cross section of the Drell-Yan process, while dσ EWIMP /dm is an EWIMP contribution to the cross section. Thus, with introducing the parameter µ, we denote the total cross section as Obviously, µ = 0 corresponds to the SM, while µ = 1 corresponds to the SM+EWIMP model. Hereafter, we use which parameterizes the amount of the correction. Note that this ratio remains unchanged even if we take into account the next-to-leading order (NLO) QCD effect because the EWIMPs affect the cross sections only through the vacuum polarization. #3 In Fig. 1, we show δ uu σ (left) and δ dd σ (right) as a function of m . The red, blue, orange, and gray lines correspond to Higgsino, wino, 5-plet fermion, and 7-plet scalar, respectively, with their masses fixed to be 1 TeV. Note that the lines sharply bend at m = 2m, which is clearer for the fermionic case. These bends are due to the fact that the pair production channel of the EWIMPs opens just above m = 2m. It is this characteristic shape that plays an important role for the discrimination of the EWIMP signal from various systematic errors (see Sec. 3.2) and also for the determination of the EWIMP mass (see Sec. 3.4). #3 When the NLO QCD effect is included, one of the initial partons can be gluon with the real emission of one jet in the final state. However, we can easily see that δ ug σ = δ uu σ and so on.  3 Analysis

Event generation
In our analysis, we first generate events for the Drell-Yan process within the SM. We assume the center of mass energy √ s = 100 TeV and the integrated luminosity L = 30 ab −1 .
For the event generation, we use MadGraph5_aMC@NLO [50,51] for the hard process, followed by parton shower and hadronization with Pythia8 [52], and detector simulation with Delphes3 [53]. We use NNPDF2.3QED with α s (M Z ) = 0.118 [54] as an input set of PDFs. As renormalization and factorization scales, we take the geometric mean of lepton transverse momenta (which we call Q). We take into account the NLO QCD effect by the [QCD] option of MadGraph5_aMC@NLO since it enhances the cross section roughly by a factor of 2 compared to the LO calculation. #4 Events with = e, µ are generated within the range of 500 GeV < m < 15000 GeV and divided into 145 bins with 100 GeV width each. We use the detector card for the FCC-hh detector installed in Delphes3 by default. The main detector effects are i) reducing the total number of events by about 20 % because of the lepton miss-identification and ii) smearing lepton momenta by a unit of O(1 − 10) GeV depending on the original momentum. The EWIMP effect can be incorporated by using δ ab σ defined in Eq. (13) and rescaling the SM event. The event number in the SM + EWIMP model in the bin m min < m < m max is theoretically estimated as where the sum runs over all events whose lepton invariant mass m obs observed by the detector simulation falls into the bin. #5 Note that the true invariant mass m true (and a, b), which is extracted from the hard process information, is used for calculation of δ ab σ (m true ) since m true = m obs due to the detector effect in general.

Statistical treatment
In order to obtain the experimental detection reach, we use the profile likelihood method. With the cross section given in Eq. (12), we calculate the theoretically expected number of eventsx i (µ) in the i-th bin by substituting δ ab σ → µδ ab σ in Eq. (14). Up to now, we have evaluated the expected number of events without taking into account systematic errors. However, in the actual experiment, the number of events may be affected, for example, by the uncertainties of the integrated luminosity and the beam energy, choices of the renormalization scale, the factorization scale, and PDF, detector efficiency and acceptance, and so on. In order to take account of the systematic errors in association with these uncertainties, we introduce nuisance parameters θ = {θ α }, by which some (or all) of the uncertainties might be absorbed. We express our theoretical prediction of the number of events in the i-th bin including the systematic errors as where f sys,i (θ) is some function of the nuisance parameters that has a property f sys,i (0) = 1.
The explicit form of f sys,i used in our analysis will be given below.
In order to discuss the detection reach, we define the following test statistic; once the number of events are experimentally measured, we introduce [55] q 0 = −2 ln L(x;θ, µ = 0) where x = {x i } with x i being the observed number of events in the i-th bin andθ and {θ,μ} are the best fit values of θ (and µ) which maximize L(x; θ, µ = 0) and L(x; θ, µ), #5 In the actual analysis, sufficient number of events are generated and each event is rescaled by some weight w in order to adjust the total event number to the corresponding value under 30 ab −1 . In this case, w should be multiplied to the right-hand side of Eq. (14).
respectively. The likelihood function is defined as Here, L (θ; σ) is the likelihood function for the nuisance parameters with σ (which we call "variance") parameterizing the possible sizes of θ in the SM. The limit σ → 0 corresponds to an optimistic case without any systematic errors. We evaluate q 0 by taking x = {x i (µ = 1)} to estimate the detection reach of the EWIMPs in the following. (For such a choice,θ = 0 andμ = 1.) We identify √ q 0 1.96 and √ q 0 = 5 as the reach at 95 % C.L. and 5σ levels, respectively. If we instead study the exclusion prospect, we may take x = {x i (µ = 0)}. Because the difference between the number of events for µ = 0 and µ = 1 is numerically small in the present case, the detection reach we will show in the following can also be regarded as the exclusion reach of EWIMPs with good accuracy (see Eqs. (12) and (14) of Ref. [55]). Now we move on to the details of our treatment of the systematic errors. We do not know the best choice of f sys,i (θ) before the start of the experiment. In this paper, we adopt what is known to work for the LHC data [56] and take the functional form of f sys,i (θ) as where p i = 2m ,i / √ s with m ,i being the median of the i-th bin lepton invariant mass. We will check that the effects of the systematic errors of our concern can be absorbed in θ with this form of f sys,i (θ). In order to estimate the variances of the nuisance parameters σ = {σ α }, we first perform the fit of the results to potential sources of systematic errors using Eq. (20). We prepare two data sets (i.e., number of events in each bin) y = {y i } and y = {y i }, both of which are obtained assuming the SM only hypothesis. Here, y is given by the original set up, while y is calculated with varying one of the followings: Then, we determine the best fit values of θ by minimizing the following χ 2 variable: Sources of systematic errors  Table 1: Standard deviation of the parameters θ derived through the fit of various systematic errors. In the last line, the standard deviations are combined for each parameter assuming that all the errors are independent with each other.
The best fit values are added in quadrature to estimate the variances. The effect of the luminosity uncertainty is estimated by varying it by ±5 %. The effects of renormalization and factorization scales uncertainties are studied by varying them to 2Q and to Q/2. The effect of the PDF choice is studied by the "systematics" module that is a built-in program of MadGraph5_aMC@NLO. We have checked that, for the sources of the systematic errors listed above, y is well fitted by our choice of f sys,i (θ). In Tab. 1, we show the variances of θ for each source of the systematic error. Assuming that all sources of the systematic errors are uncorrelated, we combine all the systematic errors; the result is shown in the last line denoted as "Total", which is used as σ in Eq. (19) when we include the systematic errors. Several comments on other possible sources of systematic errors are in order. As for the beam energy error, we could not generate events at NLO due to the lack of sufficient computational power. Instead, we checked at LO that the corresponding values of σ (assuming that the uncertainty of the beam energy is 1 %) are small enough, and hence we simply ignored it. Two of the remaining sources are the pile-up effect and the underlying event, but they may be thought of as negligible since we are focusing on the very clean signal of two energetic leptons. Another one is the effect of background processes which is not considered in our analysis. It is in principle possible to estimate its effect through the simulation and improve the analysis but here we just leave it as a future task. Yet another source is the error in the simulation of detector effect which can not be treated in our procedure. Related to this, we note here that a smooth change of the event number in general, possibly including the uncertainty related to the detector effect, could be absorbed by a minimization procedure using some fit function like in Eq. (20). On the other hand, as we will discuss below, the EWIMP signal can not be fully absorbed by the fit because of the sharp bend we mentioned before.

Detection reach
Now we show the detection reach of EWIMPs at future 100 TeV hadron colliders. In Fig. 2, we plot the value of √ q 0 as a function of the EWIMP mass. The solid and dashed lines correspond to the optimistic detection reach without systematic errors (i.e. σ → 0) and the  reach obtained using the five parameters fit, respectively. The red and blue lines show the result with Higgsino and wino in Fig. 2a, while the orange and gray lines are the result with 5-plet fermion and 7-plet scalar in Fig. 2b. Let us start with the detection reach of Higgsino and wino in Fig. 2a. We can see that the 5σ discovery reach is about 1.2 TeV and 1.6 TeV, while the 95 % lines are 2.0 TeV and 2.7 TeV for Higgsino and wino in the case without systematic errors. Taking into account the systematic errors with σ given in Tab. 1, these reaches are lowered to 0.7 TeV and 1.0 TeV for 5σ and 1.3 TeV and 1.8 TeV at 95 % C.L. The experimental reach is better for wino than Higgsino due to its large SU (2) L charge. However, the wino DM could be as heavy as 2.9 TeV [57][58][59] within the thermal freeze-out mechanism, which is beyond our sensitivity. On the other hand, for the Higgsino case, the mass of 1.1 TeV [14] favored by the relic abundance is well within our reach at 95 % C.L. even if we take account of the systematic errors.
We compare our results with those expected in the direct search at future hadron colliders [49]. Adopting the same beam energy and the luminosity (i.e., √ s = 100 TeV and L = 30 ab −1 ), the Higgsino mass reaches at 95 % C.L. are estimated as mH < 0.9 -1.4 TeV by the mono-jet search (with 2 % -1 % systematic errors) and mH < 1.1 -1.5 TeV by the disappearing track search (with 500 % -20 % uncertainty in background estimation). #6 As we have already mentioned, the disappearing track search is relevant only when the mass difference of Higgsino is small enough. In any case, the indirect search based on the precision #6 Note that the result of the disappearing track search described in [49] assumes a detector design similar to the ATLAS tracking system for the Run-2 13 TeV LHC with Insertable B-Layer [60]. measurements gives comparable or better sensitivity to Higgsino. Also, we stress that the five parameters fit used here is just an example and that we may be able to discover Higgsino around 1 TeV by reducing the systematic errors. The situation is different for wino. The mass difference of the charged and the neutral components of wino tends to be very small (∆m 165 MeV [61]). #7 Accordingly, the disappearing track search provides by far the most efficient way for the wino search with the sensitivity up to 6 TeV [49]. Yet, of course, our analysis could provide some independent information also for the wino search.
We also comment on the detection reach of the MDM scenario shown in Fig. 2b. The 5σ reaches are 3.9 TeV and 2.2 TeV for 5-plet fermion and 7-plet scalar, while the 95 % reaches are 6.2 TeV and 3.6 TeV. They are lowered to 2.6 TeV and 1.2 TeV (5σ) and 3.9 TeV and 2.1 TeV (95 % C.L.) when the systematic errors are included. If we assume the vanilla thermal freeze-out scenario, the mass should be 10 TeV for 5-plet fermion and 25 TeV for 7-plet scalar [14]. Thus, our method probes only a part of the allowed mass range for these multiplets.
Next, in order to investigate how it is important to reduce the systematic errors, we repeat the same procedure assuming that the dominant sources of the systematic errors are well understood and σ are all reduced by a factor of ten. We obtain 5σ reach as 0.8 TeV and 1.1 TeV for Higgsino and wino, both of which are improved by 100 GeV. Such an improvement is expected if, in this example, we can reduce the errors from the renormalization scale choice through the calculation of higher order effects, and from those related to PDF by obtaining more profound knowledge about PDF itself. If we can reduce the systematic errors further, the detection reach ultimately approaches to the optimistic ones.
Finally, in order to take a closer look at the significance of the peak structure (Fig. 1), we plot in Fig. 3 the contribution of each bin to the value of q 0 . The color and line style conven-#7 In terms of the low energy effective theory, at least a dimension-seven operator is required to generate the mass difference for wino other than the electroweak correction [62].
tions are the same as those in the previous figure. We can see that the most contributions come from the bins around the peak at m = 2m. This feature is clearer for the fitting based approach, where all the smooth parts of the correction are absorbed into the fit parameters, thus there is almost no contribution to q 0 from the bins other than m ∼ 2m. Note also that, for the optimistic bound, there are more contributions from the bins with lower m than those with higher m , though sometimes the cross section correction is much larger in the latter bins. This is just because of the difference of number of events in each bin, that is O(10 7 ) for 500 GeV < m < 600 GeV, while O(10 3 ) for 4900 GeV < m < 5000 GeV in our set up, for instance.

Mass determination
In this subsection, we briefly discuss the mass determination of the EWIMP after its discover, which is also possible because of the characteristic shape of the signal. Here, we assume that the underlying model is the SM + EWIMP model with µ = 1, and calculate the theoretical expectation of the number of events in each binx i as a function of the EWIMP mass m, i.e.x i =x i (m). We further include the effects of the systematic errors using x i (θ, m) ≡x i (m)f sys,i (θ). Since our purpose here is to discuss the mass determination after the discovery, we assume that the observed number of events is given by x = {x i (m true )} in the following. Then, the accuracy of the determination of the EWIMP mass is estimated by the following test statistic: whereθ and {θ,m} are the best fit values of θ (and m) which maximize L(x; θ, m fit ) and L(x; θ, m), respectively. The likelihood function is and L (θ, σ) given by Eq. (19) with σ being 0 for the optimistic case, and being the same as those in Sec. 3.2 for the case with the systematic errors. In Fig. 4, we show contours of the test statistic for Higgsino in the m true vs. m fit plane, from which we can discuss the precision of the mass determination. The purple and golden lines correspond to √ q m fit = 1 and 2, respectively. The dotted and solid lines are the results with and without the systematic errors. As shown in the figure, for light enough Higgsino, the Higgsino mass can be determined with the precision of ±(1 -100) GeV at 1σ depending on the value of the mass itself. On the other hand, the upper bound on the Higgsino mass becomes divergent when m true is equal to the discovery reach of the Higgsino at the Higgsino mass determination 1σ 2σ 1σ (sys) 2σ (sys) Figure 4: Contours of the test statistic for Higgsino in the m true vs. m fit plane. The purple and golden lines correspond to √ q m fit = 1 and 2, respectively. The dotted and solid lines are with and without the systematic errors. For Higgsino within the 2σ discovery reach, these contours can be interpreted as 1σ and 2σ errors for the mass.
corresponding C.L. (see Fig. 2), as expected from Eqs. (16) and (23). #8 (Notice that, as m → ∞, δ ab σ (m ) → 0.) Before closing this section, we comment here that the analysis proposed here can be also used to discriminate the quantum number of EWIMP. For instance, we have checked that the signal due to the existence of 1 TeV Higgsino cannot be well-fitted by wino, 5-plet fermion, nor 7-plet scalar irrespective of their masses.

Conclusion
We have studied the prospect of the indirect search of EWIMPs utilizing the Drell-Yan process at future 100 TeV hadron colliders. First, we performed an optimistic analysis neglecting all the systematic errors, and obtained the reach up to 2.0 TeV and 2.7 TeV at 95 % C.L. for Higgsino and wino, respectively. Next, we have also taken into account systematic errors by using a five parameters function fit and seen that the limits are lowered to 1.3 TeV and 1.8 TeV. These bounds for Higgsino are comparable to or better than the other search strategies that make use of direct production of Higgsino at the colliders. For wino, the analysis proposed here can also provide some information independent from the direct search. We have also shown the 95 % C.L. reach for 5-plet fermion and 7-plet scalar: 6.2 TeV and 3.6 TeV for the optimistic analysis and 3.9 TeV and 2.1 TeV for the analysis with a fitting procedure. The above values can also be interpreted as the upper limit on the EWIMP #8 For the case with systematic errors, the contours giving the upper bounds on the mass enter the regions above the discovery reach. This is due to the fact that, with systematic errors, the value of q m fit is mostly determined by the bins around m true and m fit . masses if we instead consider the exclusion prospect.
We have also applied our analysis to the determination of the EWIMP mass after its discovery. As an example, we have shown the result of the Higgsino mass determination. The Higgsino mass can be reconstructed with an uncertainty of ±(1 -10) GeV at the 1σ level, depending on the mass itself.