Probing the Dark Sector through Mono-Z Boson Leptonic Decays

Collider search for dark matter production has been performed over the years based on high pT standard model signatures balanced by large missing transverse energy. The mono-Z boson production with leptonic decay has a clean signature with the advantage that the decaying electrons and muons can be precisely measured. This signature not only enables reconstruction of the Z boson rest frame, but also makes possible recovery of the underlying production dynamics through the decaying lepton angular distribution. In this work, we exploit full information carried by the leptonic Z boson decays to set limits on coupling strength parameters of the dark sector. We study simplified dark sector models with scalar, vector, and tensor mediators and observe among them different signatures in the distribution of angular coefficients.Specifically, we show that angular coefficients can be used to distinguish different scenarios of the spin-0 and spin-1 models, including the ones with parity-odd and charge conjugation parity-odd operators. To maximize the statistical power, we perform a matrix element method study with a dynamic construction of event likelihood function. We parametrize the test statistic such that sensitivity from the matrix element is quantified through a term measuring the shape difference. Our results show that the shape differences provide significant improvements in the limits, especially for the scalar mediator models. We also present an example application of a matrix-element-kinematic-discriminator, an easier approach that is applicable for experimental data.


Introduction
The existence of dark matter (DM) is now well established. Current measurement gives a cold DM density of 25.8%, which is much significant than the 4.84% baryon density [1,2]. Despite being an essential constituent of the universe, intrinsic properties of the DM, like mass, spin and nongravitational interaction between the standard model (SM) particles are still elusive at present. Assuming that DM is weakly interacting with the SM particles, the DM annihilation cross section will be constraint by the precisely measured relic DM abundance and a weak-scale DM candidate is usually expected for consistency [3]. The WIMP DM candidate can be produced at the LHC, and its missing from detection typically leads to large missing transverse energy, resulting in mono-X signatures, where X may denote a jet [4][5][6], especially t-/b-jet [7, 8], a photon [9], a Z boson [10-12], a W boson [13,14] or a Higgs boson [15,16]. Numerous efforts have been performed at the LHC searching for the DM, many results from 13 TeV collisions are now available [5-9, 11, 12, 14-20], with strategies and benchmark models described in Ref. [21].
In this analysis, we explore the effectiveness of the Z boson leptonic decay with mono-Z signature in probing properties of the dark sector. Compared with other search channels, this channel has a relatively lower cross section and may not be the most powerful one at the stage of searching. However, precisely measured electrons and muons provide a clean signature and can be used to increase the signal feasibility. Phenomenology of this channel has been explored in Ref. [22][23][24][25][26], including higher-order QCD predictions, multivariate analysis, a search for extra dimension and effects on electron-positron colliders. LHC measurements are also available, and limits have been set on several dark sector models [11,12,27]. To better exploit the powerfulness of the lepton angular distribution, we study systematically information carried by the angular distribution and how they are affected by the dark sector.
The modeling of the dark sector can be implemented in many models. As there is no strong support for the correctness of a specific model, it is now popular to set limits on parameters of effective or simplified theories [21,[28][29][30]. Despite the simplicity, these models may not be realistic if we are not applying them in a suitable case. Either oversimplification nor overdress of the theory can lead to ineffectual results. For example, going to very high energy can result in the violation of unitarity in effective theories [21,31]. On the other hand, some features are general among models and can have less dependence on the variations of model parameters, e.g., spin and mass of the dark mediator, parity or charge conjugation parity (CP) of the couplings. If applying carefully, those effective or simplified models can help us better understand the phenomenology of the dark sector.
Motivated by this, we look for specific variables that can have discrimination power on general features of the dark sector. We consider the associated production of a Z boson and a dark mediator, where the Z boson decays to a pair of electrons or muons and the dark mediator decays to a pair of dark matter. As the dark matter is unmeasurable, the typical feature of the event is a single leptonically decaying Z boson, with p T balanced by the missing transverse momentum vector. With precisely measured electron or muon momenta, one can reconstruct the Z boson rest frame and study in detail information carried by the Z boson spin density matrix. We consider simplified models for spin-0, spin-1, and spin-2 mediators [22,[32][33][34][35]. In each case, only a few benchmark scenarios are considered with representative parameter values. For the spin-0 model, we assume the dark mediator can only weakly interact with bosons through a set of dimension-5 operators as described in Ref. [22]. In this case, the mono-Z boson channel is advantageous as a triple boson coupling is necessary for the production. If introducing couplings to the SM fermions assuming minimal flavor violation, their effects are suppressed due to proportionalities to the Yukawa couplings [33,36,37]. The spin-1 mediator model is chosen to be consistent with the one adopted in the LHC experiment [21]. A spin-2 mediator model described in Ref. [35] is also tested.
To maximally exploit the statistical power of the data, we present a framework to use the matrix element method (MEM) with a dynamical construction of event likelihood function and set unbinned limits on parameters of the dark sector [38][39][40][41][42]. We parametrize the test statistic in a way such that the sensitivity of MEM can be quantified through a term proportional to the KL-divergence of two probability density functions [43]. Limits on the coupling strengths of the dark sector models are set at 95% confidence level (CL) based on the asymptotic approximation. As the spin-2 scenarios are found to have similar angular coefficients to the one of a spin-independent spin-1 model, they are not considered in the limit setting. An example application of a matrix-element-kinematic-discriminator is also demonstrated with simulated events. This paper is organized as follows: Section 2 introduces the parametrization of lepton angular distribution. Section 3 describes computational details and presents numerical results of angular coefficients in the Collins-Soper frame. Section 4 explained the statistical method for setting limits and present results on the coupling strengths of dark sector models. Section 5 summarizes our major findings and outlooks aspects of the study.

Parametrization of lepton angular distribution
A probability density function (pdf) for a single event can be defined through the matrix element as [44] where Φ represents the Lorentz invariant phase space, in our case, a four body version Φ 4 (k l , kl, k χ , kχ) with l = e, µ and χ for the DM particle. f a (x, µ F ) corresponds to the parton distribution function of parton a, with an energy fraction of x and a factorization scale µ F . λ stands for a set of parameters of interest. The visible part of the phase space is determined through observables, while the invisible part is integrated over. The general cross section formula is written as: For the same process, it follows that the ρ(p vis |λ) is indeed a probability density function for the visible kinematics: To calculate the production of a Z boson in association with a DM mediator, we parametrize the four-momenta as follows: It is common to study the decaying lepton angular distribution in the Collins-Soper (CS) frame [45]. The Collins-Soper frame, as shown in Fig.1, is a Z boson rest frame, with the z-axis lying in a way bisects the opening angle θ ab between the beam and negative target momenta directions. In this frame, momenta of the two incoming partons become: where the x 1,2 and y Z dependences have been factorized out. Determined by these two momenta, the z-axis of this frame treats the in-and out-partons equally and tan θ ab is invariant under the longitudinal boost. This feature makes it suitable for the study of effects at finite |q T |. To avoid possible dilutions by the initial states swapped processes, we performed a rotation of π around the x-axis for events with y Z < 0 [46,47]. This rotation makes all angular coefficients distribute symmetric in y Z . In experiment, only the two decaying lepton pair are measurable, giving a set of visible variables y Z , q T , s Z , cos θ CS , φ CS , where the latter two denote polar and azimuthal angles of the charged lepton in the CS frame. We parametrize the Lorentz invariant phase space in a way such that the invisible part s Y , y Y , cos θ χ , φ χ can be integrated over: Then we factorize the decay angular distribution in terms of nine harmonic polynomials and eight angular coefficients A i , i = 0, ..., 7 [46,47]: where the polar and azimuthal angles θ, φ are measured in the CS frame. Coefficients A 5 − A 7 are parity-odd and do not contribute at tree level and are found to be very small for a Z boson production [46,47]. Therefore in this analysis, we consider only A 0 − A 4 .

Numerical results of angular coefficients in the Collins-Soper frame
As we are not directly searching for a resonance, the s Z is expected to give no sensitivity and a narrow width approximation (NWA) is applied for convenience. Apart from that, we have four observables from the Z boson decay: y Z , q T , cos θ CS , φ CS . To study the features of this four-dimensional data, we calculate angular coefficients in the y Z − q T plane for both the major background process ZZ → 2l2ν production and different dark sector models. The angular coefficients can be extracted using the method of moments [48]. In the experiment, it is more straightforward to extract from a likelihood fit [46,47].
Applying NWA for the Z boson, the cross section can be calculated through spin density matrices of the Z boson production (ρ P ) and decay (ρ D ): The Z boson production density matrix is defined in a specific range (R) of y Z − q T as follows: where ext means sum over spins and colors of all external particles other than the Z boson and averaged for the initial state ones. The decay density matrix is obtained using the Z boson decay amplitudes and parametrized similar as in Ref. [49]. The production and decay density matrices are both normalized such that the trace is one.
To obtain the amplitudes, we start from the FeynRules models implimented by authors of Ref. [22,[32][33][34][35] and use ALOHA in the MadGraph framework to generate HELAS subroutines for the helicity amplitudes [50][51][52][53][54]. In the CS frame, we choose the z-axis as spin quantization axis, hence a rotation is necessary to bring the helicity frame results to the CS frame ones. We choose the y-axis to be common for the two frames and find the opening angle between the two frames ω can be obtained through where τ Z ≡ s Z s and ω ∈ [0, π). The density matrices are then rotated according to Wigner's d-functions: where we have used the following notations: The phase space is prepared analytically, and integration is performed using BASES [55] and GNU Scientific Library. We mapped the phase space variables to increase integration efficiencies. Specifically, for a massive propagator with mass m and width Γ, the invariant mass is generated with s = m 2 + mΓ tan(x(y max − y min ) + y min ), where (3.5) and x is a uniformly generated random number. The simulation considers sin θ W = 0.23129, m Z = 91.1876 GeV, Γ Z = 2.4952 GeV and α(m Z ) −1 = 127.95 [1]. The W boson mass is obtained through m Z cos θ W , assuming ρ parameter equals to one. The α S is chosen to be consistent with the one in the parton distribution functions (PDF). We use PDF set NNPDF23 [56] with α S (m Z ) = 0.130 at leading order. The factorization scale is set to be equal to the Z boson transverse energy Cross sections in this section consider the visible Z boson decays to electrons and muons with NWA and Br(Z → l + l − ) = 6.73% [1]. The advantage of our program is that high statistical accuracy can be achieved through a direct integration. To validate our program, we checked our angular coefficients through toy measurements based on MadGraph5 aMC@NLO (MG5) generated events.

SM ZZ → 2l2ν background
The SM ZZ → 2l2ν production is the major background of our DM search. It has a similar final state signature as the signal process, as depicted in Fig. 2. Hence we first take a look at the Fig. 3 for the angular coefficients of this process. In general, the angular coefficient A 0 measures the difference between longitudinal and transverse polarizations, and it looks more longitudinal at high q T . The coefficient A 4 measures forward-backward asymmetry, the Z boson looks more like left-handed in the forward region. The A 2 measures the interference between the transverse amplitudes and the A 1,3 measures the interference between transverse and longitudinal. Figure 2. Representative Feynman diagrams of the SM ZZ → 2l2ν production.   We consider a simplified model with a scalar s-channel mediator as described in Ref. [22]. The dark sector model is constructed as follows:

Spin-0 mediator
whereṼ µν = 1 2 ǫ µνρσ V ρσ is the dual field strength tensor of V field, Λ is a high energy scale. As discussed in Ref. [22], this operator can be induced by a fermion loop graph with heavy fermion integrated out. Signature of this model is very different from the SM ZZ → 2l2ν process, the dark mediator is emitted from the SM gauge bosons as depicted in Fig. 4. We consider three benchmark scenarios of the parameters labeled by S0 a,b,c . As our angular distributions are more sensitive to changes in couplings, we fix the mass of dark matter m χ = 10 GeV and the mass of the mediator m Y 0 = 1000 GeV. The angular distributions won't be changed drastically as long as 2m χ is much smaller than m Y 0 . The parameter values and inclusive cross sections are listed in Table 1.
Angular coefficients of the benchmark scenatios S0 a,b,c are shown in Fig. 5, Fig. 6 and Fig.7 respectively. Comparing to the SM ZZ → 2l2ν, the dark matter signal is produced with much higher q T and have very different angular coefficients distributions, e.g., more transverse at low q T . The S0 a and S0 b can be distinguished from A 0 , A 2 , where the y Z dependences are very different. In the case of S0 c , Y 0 couples to weak bosons like a Higgs boson and cannot perturb the coupling structure with the Z boson production. Consequently, the A 0 , A 1 and A 3 in the CS frame are all zero hence are not shown in the figure. Benchmark

Spin-1 mediator
We consider the same dark sector with a spin-1 mediator as in the LHC experiment [21] with the following interactions of the dark sector: The masses of the dark matter and the mediator are chosen to be the same as in the spin-0 model. A sound discussion of the impact of the choice of masses is available in the Ref. [21]. Since our analysis is more suitable for testing couplings, we consider benchmark scenarios as listed in Table 2. The signal signature is close to the SM ZZ → 2l2ν process, as shown in Fig. 8, and we include here the SM ZZ → 2l2ν as a special case with zero coupling for comparison. The S1 b and S1 c project out the right-and left-handed part the Z-q-q couplings. Since the magnitude of the left-handed couplings are larger than the one of the right-handed, cross section of the S1 c scenario is found to be much larger than the S1 b scenario.
Angular coefficients of the benchmark scenarios S1 a,b,c are shown in Fig. 9, Fig. 10 and Fig. 11 respectively. Comparing with the SM ZZ → 2l2ν and spin-0 dark sector models, A 0 of the spin-1 models are found to be very significant. Among the three scenarios, most signatures look similar, but A 3 and A 4 take different signs between the S1 b and S1 c . Hence the A 3 and A 4 can be used to quantify the parity violation of the dark sector.
Benchmark S1 a S1 b S1 c S1 0 Spin independent Right handed Left handed SM (ZZ → 2l2ν) 56.3 55.9 55.9 -Cross section (fb) 2.50 0.533 4.50 239 Table 2. Benchmark scenarios with a spin-1 mediator.  Figure 9. Angular coefficients A 0 − A 4 and the y Z − q T differential cross section of the benchmark scenario S1 a .  Figure 11. Angular coefficients A 0 − A 4 and the y Z − q T differential cross section of the benchmark scenario S1 c .

Spin-2 mediator
The dark sector with a spin-2 mediator is also tested. We consider a model as described in the Ref. [35], with benchmark scenarios listed in Table 3. The masses of the dark matter and the mediator are also chosen to be the same as in the spin-0 model. Despite an increase of complexity in the computation, we found the angular coefficients look similar to the benchmark scenario S1 a . We show only the angular coefficients of the benchmark scenario S2 a in Fig. 12. Some visible differences from the S1 a can be observed from the A 0 and A 2 distributions. Since we do not measure the DM, the angular coefficients of S2 b,c are found to be very close to the ones of S2 a .  Table 3. Benchmark scenarios with a spin-2 mediator. Angular coefficients of the three scenarios look all the same.

Setting limits on the coupling strength parameters of dark sector models
In Section 3, we have shown that angular coefficients of the benchmark dark sector models can have distinct signatures from the SM ZZ → 2l2ν background process in the y Z − q T plane. In this section, we take advantage of these signatures and set limit on the coupling strength parameter λ of each dark sector model, based on observables x = (y Z , q T , cos θ CS , φ CS ). The invisible part (y Y , s Y , cos θ χ , φ χ ) was integrated out to construct pdfs, as described in Section 2.

Statistical method
With the pdfs of the signal and background processes obtained through MEM, one can construct an unbinned likelihood function over N events in the data sample [57]: where ρ s (x, λ) and ρ b (x) represent pdfs of the signal and background, S(λ, θ) and B(θ) corresponding to the expected signal and background yields. The θ represents the full set of nuisance parameters with pdf ρ(θ), which are designed to incorporate systematic uncertainties.
To set limits on the parameters λ, we compare the compatibility of the data with the λ fixed and λ floated hypotheses and construct a test statistic based on the profile likelihood ratio: According to the Wilk's theorem, this test statistic satisfies the χ 2 distribution of the same degrees of freedom as λ in the large sample limit [58]. One can, therefore, set limits on the λ through a parameter space scan and cut on the −2 ln ∆L values.
Neglecting pdf of the nuisance parameters, it follows that For setting limits on λ, we assume that there is a single dataset in agreement with λ = 0. In the large sample limit, we have: where the first term is a test statistic for simple counting experiment and the second term is proportional to N and a KL-divergence [43]. As the KL-divergence measures the difference of the pdfs ρ(x|λ) and ρ(x|λ = 0), it quantifies the powerfulness of the MEM. For simplicity, we will call the first term as normalization term and the second one as KL-divergence term.
In our study, the likelihood function is prepared by BASES numerical integration with HELAS subroutines for the helicity amplitudes. The evaluation of the KL-divergence term is performed using a plain integration provided by the GNU Scientific Library. We validate our program by checking the normalizations of all the constructed pdfs and by comparing the angular coefficients and cross sections of all involved processes with the MG5. See more information in Appendix A.

Background modeling and event selections
To make our limits more realistic, we consider a few selections -marked as BL selectionsas listed in Table 4 to capture major detector acceptance effects for the processes involved. The values of these selections are set refering to recent 13 TeV LHC measurements [11,12]. There are several additional selections considered in experiments to improve the signal feasibility, e.g., jet counting, 3 rd -lepton veto, top quark veto, and ∆φ ll,p miss . These selections reject most background from misidentification but lead to different acceptance efficiencies for different processes. Without detector simulation, we determine the event rate according to the CMS results (Table 3 of Ref. [11]), with an ancillary A · ǫ incorporating the additional selections in the experiment and a scale factor normalizing to 150 fb −1 data. The signal dark matter processes are assumed to have the same ancillary A · ǫ as the SM ZZ→ 2l2ν process. Table 4. Selections considered in our computations (BL-selections), where l = e, µ. Additional selection requirements are considered in experiments to improve the signal feasibility. Their effects are included through an Ancillary A · ǫ.

Variable Requirements
Our background pdf is constructed based on components summarized in Table 5. Apart from the non-resonant-ll background, which is constructed using only the phase space, other components are built using matrix elements. The WZ→ 3lν matrix element assumes W→ eν, where the electron is not identified by a detector. The Z/γ * → l + l − is estimated with matrix element of the Z→ l + l − plus one jet production, phase space of this process reduces to three final state particles.
[11]. The number of events has been translated into 150 fb −1 data.
In the presence of selections, angular coefficients can be distorted. Fig. 13 shows the angular coefficients A 0 − A 4 for the background only hypothesis. Irregular distributions on the boundaries are mainly caused by the selections on |η l | and ∆R ll . With the coupling strength at our expected limit, the presence of signal can only perturb the shapes of the background only ones.  Figure 13. Angular coefficients A 0 − A 4 in the CS frame and y Z − q T differential cross section for background only hypothesis. Selections in Table 4 have been applied and cause irregular shapes in kinematic boundaries.

Limits on the coupling strength parameters of the dark sector models
In our dark sector models, it is necessary to have two couplings: one for the interaction with SM particles, one for the DM decay. For conciseness, we assume that both couplings in the benchmark model are scaled by a strength parameter λ. This assumption makes the cross sections change with two orders severer in couplings than ones for limits of a single coupling. We compare the upper limits set from the normalization term −2 ln Poisson and from the KL-divergence term 2N · D(ρ(x|0)||ρ(x|λ)) in Fig. 14 for the S0 benchmark scenarios and in Fig. 15 for the S1 benchmark scenarios. The shapes provide significant improvements in all cases. The KL-divergence terms drive the final limits for the S0 benchmark scenarios and are close to the normalization terms in the S1 benchmark scenarios.
68% CL. limit 95% CL. limit 99% CL. limit c S1 Figure 15. Upper limits on the coupling strength parameters of the S1 benchmark scenarios.
We provide in Table 6 95% CL upper limits of the strength parameters. In our evaluation, the numerical uncertainty of the normalization terms can be easily made negligible. However, the evaluation of the KL-divergence terms can be computationally expensive. It takes us roughly 700 × 6 CPU hours, functioning at about 2.4 GHz, for us to obtain 30%-50% uncertainties on the KL-divergence terms around the limit values. The signal cross sections at the limit values are also reported. Since a counting experiment calculate limits based on signal background yields, the results from the normalization term are almost the same. The ones from the KL-divergence terms, however, depend on the shape difference between the signal and background. As the KL-divergence is a measure the shape difference, a lower cross section means a larger the difference in shape. These quantitative results are in agreement with qualitative features of the angular coefficients among models provided in Section 3.

Benchmark
S0 a S0 b S0 c S1 a S1 b S1 c Limit from the normalization term (λ 1 ) 4.  Table 6. Upper limits on the coupling strength parameters of the dark sector models at 95% CL, with signal cross sections at the limit values.

Example application of MEKD
Our computation considered only parton level matrix element at leading order (LO). We comment that there are already efforts to extend the MEM to Next-to-Leading Order (NLO) [59] and incorporates parton shower effects [60]. There is an easier approach to exploit the LO matrix elements, called the matrix element kinematic discriminator (MEKD) [41,61,62]. This method construct a variable named MEKD that can be calculated for events with required observables. By construction, it utilizes the matrix element and can be used to distinguish the signal and background. The advantage of this method is that detector effects and theoretical uncertainties in the construction of likelihood function is independent of the application. Based on the pdfs defined as in Eq. 2.1 of the signal and combined background, we define the MEKD as: where x = (y Z , q T , cos θ CS , φ CS ) and the invisible part has been integrated out. Then we use the MG5 program to generate events for the applications. For the LO simulations, we consider the same setup as has been used in our program. For the NLO simulations, we consider NNPDF23 nlo with default renormalization and factorization scales, defined as the sum of the transverse masses divided by two of all final state particles and partons. Negatively weighted events in the NLO simulations have been incorporated consistently.
The Fig. 16 stacks MEKD distributions of both signal and backgrounds. On the left plot, all of the processes are generated with LO accuracy (NLO in QCD for Z(→ l + l − )+jet). The signal considers S0 a benchmark model with λ = 3.5. We multiplied the signal yield by a factor of five for a better demonstration. The Non-resonant-ll process is expected to be obtained from data-driven in the experiment. We mimic its contribution by using a tt(→ 2l2ν2b − jets) sample. The right plot replaces the SM ZZ→ 2l2ν, WZ(→ eν2l) and Z(→ l + l − )+jet with simulated events at NLO accuracy. In both cases, the MEKD shows very nice discrimination power on the signal and background. It is made clear that NLO simulated events are applicable, with a reasonable loss of sensitivity.

Summary
In this paper, we have exploited the Z boson leptonic decay information to probe the dark sector with a scalar, vector, and tensor mediators. We obtained angular coefficients of the SM ZZ → 2l2ν background and benchmark scenarios of the dark sector models in the y Z −q T plane. Our results show that the angular coefficients A 0 −A 4 behave very differently between the SM ZZ → 2l2ν process and the dark sector signal processes. The angular coefficients among dark sector models of spin-0 and spin-1 mediators are also found to be different from scenario to scenario. Specifically, the angular coefficients have sensitivities on the parity violation of the spin-1 model and the CP-violation of the spin-0 model. The Events / bin 10 2 10 3 10 5 Events / bin 10 2 10 3 10 5 Figure 16. Example MEKD distributions with MG5 generated events. The left plot is obtained with simulated events at LO accuracy. The right plot considers events of Z(→ l + l − )+jet processes at NLO accuracy. The signal considers S0 a benchmark model with λ = 3.5. We multiplied the signal yield by a factor of five for a better demonstration. angular coefficients in the spin-2 model are found to be similar to the spin-independent scenario of the spin-1 model but still have minor differences.
To quantify the shape information that can be used for the search of dark sectors, we consider unbinned fits to the four-dimensional y Z − q T − cos θ CS − φ CS distributions based on dynamically constructed matrix element likelihood functions and set 95% CL upper limits on the coupling strength parameters of the spin-0 and spin-1 benchmark scenarios. To be realistic, we emulate the acceptance and efficiency effects referring to the 13 TeV LHC measurement [11,12]. To make our framework concise, we obtained all the results using asymptotic approximation without event generation.
Our evaluated KL-divergence term quantifies the shape effect in each case. The obtained results demonstrate significant improvements in the limits, especially on the S0 benchmark models. For easier usage of experimental data, we provide an example application of MEKD with simulated events. We show that our MEKD constructed with LO matrix elements are applicable for NLO events and preserves good discrimination power on the signal and background. We expect this kind of MEKDs to be useful for exploiting the lepton angular distributions in experimental analyses. and No. 1157500, and by a short-term internship program from the graduate school of Peking University.

A Cross checks with the MG5 program
To make the MG5 results comparable, we implemented similar setups as described in the paper. These include coupling constants, the choice of PDF set, renormalization and factorization scales, Breit-Wigner cutoff, and BL selections as described in Table 4 of the paper. The Table 7 compares our results with the MG5 ones with one on-shell Z boson in the final states. For all the cases, the differences lie within statistical uncertainty. The Table 8 compares our results with the MG5 with the Z boson leptonicalled decayed. Our program considered all the BL-selections with NWA, while the MG5 ones replace the NWA with |m ll − m Z | < 15 × Γ Z . This replacement leads to slightly smaller MG5 cross sections comparing to ours, but in general, the differences are not large. Normalizations of the signal and all of the background pdfs are also checked to be consistent with one.