Exploring properties of long-lived particles in inelastic dark matter models at Belle II
KKIAS-P21001
Exploring properties of long-lived particles ininelastic dark matter models at Belle II
Dong Woo Kang , P. Ko , Chih-Ting Lu School of Physics, KIAS, Seoul 130-722, Republic of Korea (Dated: January 8, 2021)
Abstract
The inelastic dark matter model is one kind of popular models for the light dark matter (DM)below O (1) GeV. If the mass splitting between DM excited and ground states is small enough, theco-annihilation becomes the dominant channel for thermal relic density and the DM excited statecan be long-lived at the collider scale. We study scalar and fermion inelastic dark matter modelsfor O (1) GeV DM at Belle II with U (1) D dark gauge symmetry broken into its Z subgroup. Wefocus on dilepton displaced vertex signatures from decays of the DM excited state. With the helpof precise displaced vertex detection ability at Belle II, we can explore the DM spin, mass andmass splitting between DM excited and ground states. Especially, we show scalar and fermion DMcandidates can be discriminated and the mass and mass splitting of DM sector can be determinedwithin the percentage of deviation for some benchmark points. Furthermore, the allowed parameterspace to explain the excess of muon ( g − µ is also studied and it can be covered in our displacedvertex analysis during the early stage of Belle II experiment. a r X i v : . [ h e p - ph ] J a n . INTRODUCTION The nature of dark matter (DM) in our Universe is still a great mysterious issue. As weknow, DM plays a major role in the structure formation [1], and its abundance is about 5.5times larger than the ordinary matter in the present Universe [2]. However, even all robustevidences to support the existence of DM until now are only connected with gravitational in-teractions. Still it is believed that there are non-negligible couplings between Standard Model(SM) particles and DM in addition to gravitational interactions. Especially, the Weakly In-teracting Massive Particle (WIMP) DM candidate with the freeze-out mechanism has beenoverwhelming in both theoretical and experimental communities during the past decades asshown in Ref. [3] and references therein. This kind of WIMP DM particles can be searchedwith direct, indirect detections and also at colliders experiments. The null signal result fromdirect detection provides severe bounds on the cross section of DM and nuclear scatteringabove a few GeV [4–6]. For instance, the spin independent (spin dependent) DM and nuclearscattering cross section is restricted to σ SI (cid:46) − cm [4] ( σ SD (cid:46) × − cm [6]) for M DM = 100 GeV. Nevertheless, lighter DM candidates are still less constrained because therestrictions of fine energy threshold are required from these direct detection experiments.Therefore, those dark sector models with MeV to GeV DM candidates become more andmore popular for phenomenological studies and new experimental searches [7, 8].For the GeV scale DM searches, the high-intensity machines such as the BESIII [9],BaBar [10], Belle II [11, 12] and also the fixed target experiments [13] can be more powerfulthan the high-energy machines such as the Tevatron and LHC [14]. On the other hand, DMis not the lonely particle in the dark sector for most of DM portal models [15–18]. In order toconnect the SM sector with the dark sector, a mediator is required. Among these mediatorcandidates, the dark photon via the kinetic mixing portal is an attractive type [19] for thedark sector with Abelian gauge symmetry. In order to be consistent with the relic densitybound, it’s natural for both DM and dark photon in a similar energy scale. Therefore,the searches for them at high-intensity machines are in full swing and relevant constraintscan be set up as shown in Ref. [20, 21]. Even though the original purpose for these high-intensity machines is to study the properties of J/ψ and B mesons, we can also apply themto explore some models with light DM candidates via the mono-photon signature [20–23].Furthermore, as pointed out in Ref. [20, 21, 24], we can not only study the mono-photon2ignature, but also the displaced DM signature at B-factories for the inelastic DM models.Especially, searches for displaced DM signature can cover some parameter space which theinvisible DM signature cannot reach at high-intensity machines. Hence, we will focus onexploring the properties of long-lived particles in dark sectors with U (1) D gauge symmetryat Belle II in this work.The inelastic (or excited) DM models with extra U (1) D gauge symmetry [25–28] is oneof the most popular dark sector models with light DM candidates . There are at least twostates in the dark sectors and there is an inelastic transition between them via the new U (1) D gauge boson. If the mass splitting between these two states is small enough, the co-annihilation channel could be the dominate one of DM relic density in early Universe [29]. Itis one of the unique features of this sort of models. The co-annihilation production for lightDM via thermal freeze out is still consistent with the Cosmic Microwave Background (CMB)constraint for the amount of parameter space [28, 30]. On the other hand, the constraintfrom DM and nuclear inelastic scattering is much weaker than the elastic one in the directdetection experiments . It makes more allowed parameter space can be explored in theinelastic DM models. Besides, there are also other rich phenomenons in this model. Forexample, it is possible to explain the muon g − , / , , / In Ref. [25], the U (1) D symmetry is explicitly and softly broken by a dim–2 operator for scalar DM, andthere is no dark Higgs boson there. On the other hand, in the models considered in this paper in thatthere is dark Higgs boson which plays important roles in DM phenomenology. The bounds from direct detection experiments are weak for ∼ O (1) GeV DM. Hence we can safely ignorepossible contributions to the DM-nucleon elastic scattering from the SM Higgs boson and dark Higgsboson in the direct detection experiments.
3n Ref. [44, 45] and collider experiments in Ref. [46–60]. Especially, using the displacedvertex information for DM mass determination starts from Ref. [57]. Thanks to the precisetrack resolution in the inner detector of ATLAS and CMS experiments, the DM mass re-construction can be studied in some BSM models of the specific cascade processes involvinglong-lived particles. Furthermore, the resolutions of displaced vertex, charged lepton andphoton momentum can be even improved at detectors of Belle II experiments [11]. There-fore, our goal in this paper is to explore the spin and mass properties of inelastic DM modelswith dilepton displaced vertex signatures at Belle II. As we will see in Sec. IV, the scalarand fermion inelastic DM models can be well discriminated and the DM mass and masssplitting between DM excited and ground states can be determined within the percentageof deviation for those benchmark points.The organization of this paper is as follows. We first review scalar and fermion inelasticDM models with U (1) D gauge symmetry in Sec. II. Both analytical representations andnumerical results for cross sections of relevant signal processes are displayed in Sec. III.We also point out how to distinguish scalar and fermion inelastic DM models via the sizeand kinematic distribution of their cross sections in this section. Detailed simulations andmethods to determine DM mass, and mass splitting between DM excited and ground statesare shown in Sec. IV. Finally, we conclude our studies in Sec. V. Some supplemental formulaefor Sec. III can be found in the Appendix A and B. II. INELASTIC DARK MATTER MODELS
The scalar and fermion inelastic (or excited) DM models with U (1) D gauge symmetryare reviewed in this section. After the spontaneous symmetry breaking (SSB) of this U (1) D gauge symmetry, we expect the accidentally residual Z symmetry, φ → − φ (scalar) or χ → − χ (fermion), can be left such that φ or χ are stable and become DM candidatesin our Universe. A. The scalar model [26, 28]
We consider a dark sector with two singlet complex scalars Φ and φ = ( φ + iφ ) / √ φ are charged under4 (1) D , but neutral of the SM gauge symmetry. The U (1) D charges for them are assigned as Q D (Φ) = +2 and Q D ( φ ) = +1. Besides, the SM-like Higgs doublet and other SM particlesdo not carry U (1) D charges.The scalar part of the renormalizable and gauge invariant Lagrangian density is L scalar = | D µ H | + | D µ Φ | + | D µ φ | − V ( H, Φ , φ ) , (1)with D µ H = ( ∂ µ + i g σ a W aµ + i g (cid:48) B µ ) H, (2) D µ Φ = ( ∂ µ + ig D Q D (Φ) X µ )Φ , (3) D µ φ = ( ∂ µ + ig D Q D ( φ ) X µ ) φ, (4)where W aµ , B µ , and X µ are the gauge potentials of the SU (2) L , U (1) Y and U (1) D with gaugecouplings g , g (cid:48) and g D , respectively. The σ a is the Pauli matrix and a runs from 1 to 3. Thescalar potential in Eq.(1) is given by V ( H, Φ , φ ) = − µ H H † H + λ H ( H † H ) − µ Φ ∗ Φ + λ Φ (Φ ∗ Φ) − µ φ φ ∗ φ + λ φ ( φ ∗ φ ) + ( µ Φ φ Φ ∗ φ + H.c. )+ λ H Φ ( H † H )(Φ ∗ Φ) + λ Hφ ( H † H )( φ ∗ φ ) + λ Φ φ (Φ ∗ Φ)( φ ∗ φ ) , (5)where all parameters are assumed to be real for simplicity.We then expand H, Φ fields around the vacuum with the unitary gauge, H ( x ) = 1 √ v + h ( x ) , Φ( x ) = 1 √ v D + h D ( x )) , (6)and highlight some important relations in the below : • The four-points and off-diagonal interactions of the DM sector and new gauge bosoncan be obtained from the | D µ φ | term: | D µ φ | = ( ∂ µ φ ) + ( ∂ µ φ ) + g D X µ X µ ( φ + φ ) + g D X µ ( φ ∂ µ φ − φ ∂ µ φ ) (7) • The DM sector mass splitting and dark Higgs, DM sector trilinear interactions arederived from µ Φ φ Φ ∗ φ + H.c. terms: µ Φ φ Φ ∗ φ + H.c. = µ Φ φ ( v D + h D ( x ))( φ − φ ) . (8)5 Finally, the φ , φ masses and their mass splitting can be represented as M φ , = (cid:114)
12 ( − µ φ + λ Hφ v + λ Φ φ v D ) ∓ µ Φ φ v D , (9)and ∆ φ ≡ M φ − M φ = 2 µ Φ φ v D M φ + M φ . (10)Notice the interaction states ( h, h D ) would be rotated to the mass states ( h , h ) via themixing angle θ . The sin θ is chosen to be small enough in our analysis such that it isconsistent with the SM-like Higgs boson data at the LHC .Finally, since all SM fermions don’t carry U (1) D charges, the only way for the new X µ boson and SM fermions to interact is via the kinetic mixing between B µν and X µν . TheLagrangian density of this part can be represented as L X,gauge = − X µν X µν − sin (cid:15) B µν X µν , (11)where (cid:15) is the kinetic mixing parameter between these two U (1)s. If we apply the linearorder approximation in (cid:15) , the extra interaction terms for SM fermions and Z (cid:48) boson can bewritten as L Z (cid:48) ff = − (cid:15)ec W (cid:88) f x f f /Z (cid:48) f, (12)where c W is the weak mixing angle and x l = − x ν = 0, x q = or − depending on theelectrical charge of quark. The dark photon mass can be approximated as m Z (cid:48) (cid:39) g D Q D (Φ) v D . (13)Notice the correction from the kinetic mixing term is second order in (cid:15) which can be safelyneglected here. B. The fermion model [27, 28]
Here we consider a dark sector with the singlet complex scalar Φ and Dirac fermion χ as the dark Higgs and dark matter sectors, respectively. Both Φ and χ are charged under U (1) D , but neutral of the SM gauge symmetry. The U (1) D charges for them are assigned as The long-lived dark Higgs phenomenology study in the inelastic DM model at Belle II can be found inRef. [61]. D (Φ) = +2 and Q D ( χ ) = +1. Again, the SM-like Higgs doublet and other SM particlesdo not carry U (1) D charges.The scalar part of the renormalizable and gauge invariant Lagrangian density is L scalar = | D µ H | + | D µ Φ | − V ( H, Φ) , (14)where D µ H and D µ Φ are the same in Eq.(2) and (3). The scalar potential in Eq.(14) isgiven by V ( H, Φ) = − µ H H † H + λ H ( H † H ) − µ Φ ∗ Φ + λ Φ (Φ ∗ Φ) + λ H Φ ( H † H )(Φ ∗ Φ) , (15)where all parameters are assumed to be real for simplicity.The Lagrangian density of dark matter sector part is L χ = χ ( i / ∂ + g D / X − M χ ) χ − ( f χ c χ Φ ∗ + H.c. ) , (16)where f is assumed to be a real parameter. In order to decompose the Dirac fermion χ intoa pair of two independent Majorana fermions, χ and χ , we set χ , ( x ) = 1 √ χ ( x ) ∓ χ c ( x )) . (17)After expanding H, Φ fields around the vacuum with the unitary gauge as shown in Eq.(6),the Eq.(16) can be written as L χ = 12 χ ( i / ∂ − M χ ) χ + 12 χ ( i / ∂ − M χ ) χ − i g D χ / Xχ − χ / Xχ ) − f h D ( χ χ − χ χ ) , (18)where χ , χ masses and their mass splitting can be represented as M χ , = M χ ∓ f v D , (19)and ∆ χ ≡ ( M χ − M χ ) = 2 f v D . (20)Finally, the Z (cid:48) boson mass and its interactions with SM fermions are the same withEq.(13) and (12). 7 . The target parameter space and decay width of φ ( χ ) In this work, we focus on the scenario m Z (cid:48) > M φ + M φ ( M χ + M χ ) such that the Z (cid:48) → φ φ ( χ χ ) decay mode is kinematically allowed and becomes the dominant one. On theother hand, we are interested in the co-annihilation dominant channel for DM relic densityin early Universe, so we restrict ourselves to the compressed mass spectrum with ∆ φ,χ < . M φ ,χ in inelastic DM models [21]. Finally, we concentrate on M φ ,χ >
100 MeV suchthat the parameter space is free from the Big Bang Nucleosynthesis (BBN) constraint [62, 63] . The φ ( χ ) is the only DM candidate in our analysis. If the mass splitting ∆ φ,χ is notignorable compared with the M φ ( χ ) , φ ( χ ) is not stable and will decay to φ ( χ ) and aSM fermion pair via the off-shell Z (cid:48) . The full analytical formulas for the total width of φ ( χ ) three-body decay can be found in Appendix B of Ref. [64]. In the mass range of ourinterest, there are three kind of φ ( χ ) decay modes , φ ( χ ) → φ ( χ ) e + e − , φ ( χ ) µ + µ − and φ ( χ ) π + π − . In our numerical calculations, the total decay width of φ ( χ ) is auto-matically calculated in Madgraph5 aMC@NLO [66]. For the partial decay width of φ ( χ ) → φ ( χ ) π + π − , we rescale the partial decay width of φ ( χ ) → φ ( χ ) µ + µ − with themeasured R ( s ) values in Ref. [67]. III. THE RELEVANT CROSS SECTIONS
In this section, we show both analytical representations and numerical results for crosssections of e + e − → φ φ ( χ χ ) and e + e − → φ φ ( χ χ ) γ processes. Especially, the methodto distinguish spin-0 and spin-1/2 DM candidates in the inelastic DM models are also dis-cussed here. Notice that the DM mass can be extended to lower regions if the light dark Higgs is included as shown inRef. [28]. Since the light Z (cid:48) is isosinglet, it can mix with ω meson and decay to three pions, Z (cid:48) → π + π − π , ifkinematically allowed. However, according to Fig.1 in Ref. [65], the contributions from Z (cid:48) → π + π − π mode are important only in the adjacent regions of m Z (cid:48) ≈ m ω = 0 .
782 GeV. Hence, we don’t specificallystudy this regions in our analysis. . The analytical representations The differential cross section for e + e − → φ φ via the s-channel Z (cid:48) can be represented as dσ ( e + e − → φ φ ) dt = (cid:15) e g D πs (cid:2) ( M φ − t )( M φ − t ) − st (cid:3) [( s − m Z (cid:48) ) + m Z (cid:48) Γ Z (cid:48) ] , (21)where s , t are Mandelstam variables and Γ Z (cid:48) is the total width of Z (cid:48) boson. In order to studythe differential angular distributions of cross sections, we transfer from dσ/dt to dσ/d cos θ ,where θ is the polar angle of φ and it is defined as the direction of φ relative to the positronbeam direction.The dσ/d cos θ in the centre-of-mass (CM) frame can be simply written as dσ ( e + e − → φ φ ) d cos θ (cid:12)(cid:12)(cid:12)(cid:12) CM = 164 π (cid:15) e g D E CM [( E CM − m Z (cid:48) ) + m Z (cid:48) Γ Z (cid:48) ] ξ / (1 − cos θ ) , (22)and ξ = (cid:115) − M φ + M φ ) E CM + ( M φ − M φ ) E CM , (23)where E CM is the centre-of-mass energy. We can find Eq.(22) is always proportional to 1 − cos θ . Consequently, it indicates if the s-channel Z (cid:48) is on-shell produced, it is longitudinallypolarized (helicity = 0).However, the formula for dσ/d cos θ in the Belle II laboratory (LAB) frame is moretedious, so we show it in Eq.(A3) of the Appendix A. As we will see the numerical results,the differential angular distributions of cross sections in the LAB frame are just shifted tothe initial electron beam direction for the Belle II machine compared with the ones in theCM frame.The differential cross section for e + e − → χ χ via s-channel Z (cid:48) can be represented as dσ ( e + e − → χ χ ) dt = (cid:15) e g D πs (cid:2) s + 2( t − M χ )( t − M χ ) + s (2 t − ( M χ − M χ ) ) (cid:3) [( s − m Z (cid:48) ) + m Z (cid:48) Γ Z (cid:48) ] , (24)and dσ/d cos θ in the CM frame can be written as dσ ( e + e − → χ χ ) d cos θ (cid:12)(cid:12)(cid:12)(cid:12) CM = 132 π (cid:15) e g D E CM [( E CM − m Z (cid:48) ) + m Z (cid:48) Γ Z (cid:48) ] × (cid:20) (1 − ( M χ − M χ ) E CM + 4 M χ M χ E CM ) ξ + ξ / cos θ (cid:21) , (25)where ξ can be found in Eq.(23) with M φ , → M χ , . For the M χ , → θ , which is well known. It shows if the s-channel Z (cid:48) is on-shell9roduced, it is transversely polarized (helicity = ± ). However, the mass effects of M χ , spoil parts of this property. Similarly, we show the formula for dσ/d cos θ in the LAB framein Eq.(A4) of the Appendix A. Again, there is a skewing behavior for the differential angulardistributions of cross sections in the LAB frame compared with the ones in the CM frame.The initial state radiation (ISR) photon is an useful trigger for signals with missing energyor soft objects at lepton colliders. Especially, for the process with on-shell Z (cid:48) production,the ISR photon is used not only for background rejection, but also for Z (cid:48) invariant massreconstruction. To include the ISR photon in e + e − → φ φ ( χ χ ) process, the differentialcross section can be written as [68] dσ ( e + e − → φ φ ( χ χ ) γ ) dzd cos θ (cid:48) (cid:39) P ( z, cos θ (cid:48) ) (cid:98) σ ( e + e − → φ φ ( χ χ )) , (26)and the splitting kernal P ( z, cos θ (cid:48) ) is P ( z, cos θ (cid:48) ) = απ − z ) z θ (cid:48) , (27)where α is the fine structure constant, z = E γ E CM / is the energy fraction of ISR photonfrom the initial electron/positron and θ (cid:48) is the polar angle of ISR photon and it is definedas the direction of γ to the positron beam direction. The differential form of (cid:98) σ ( e + e − → φ φ ( χ χ )) can be found in Eq.(21) and (24). In the CM frame, if we assign the fourmomentum of the initial electron/positron and ISR photon as p e ± = ( E CM / , , , ± E CM / p ISR = ( E γ , p x,γ , p y,γ , p z,γ ), the four momentum of electron/positron after the radiationis p e (cid:48)± = ( E CM − E γ , − p x,γ , − p y,γ , ± E CM / − p z,γ ). Let’s consider two limit scenarios tointuitively catch up the behavior of d (cid:98) σ ( e + e − → φ φ ( χ χ )) /d cos θ distributons. First, ifthe ISR photon is soft ( z ≈ p e (cid:48)± ≈ p e ± and d (cid:98) σ ( e + e − → φ φ ( χ χ )) /d cos θ distributionsare close to Eq.(22) and (25). On the contrary, if the ISR photon takes almost all energyfrom the initial electron/positron ( z → φ ( χ ) is obvious in the forward or backwarddirection. Since d (cid:98) σ ( e + e − → φ φ ( χ χ )) /d cos θ distributions are highly dependent on the z parameter of the ISR photon, the e + e − → φ φ ( χ χ ) γ processes are not ideal to distinguishscalar and fermion inelastic DM models.Finally, we focus on the parameter space m Z (cid:48) > M φ ( χ ) + M φ ( χ ) in this work, so only thethree-body decay of φ ( χ ) via the off-shell Z (cid:48) is possible. The full analytical representationsof e + e − → φ φ ( χ χ ) → φ φ ( χ χ ) f f are shown in the Appendix B.10 Z ' = M ϕ χ ( Solid ) m Z ' = M ϕ χ ( Solid ) Δ ϕ χ = M ϕ χ , α D = ϵ = Δ ϕ χ = M ϕ χ , α D = ϵ = e + e - →χ χ e + e - →χ χ e + e - →ϕ ϕ e + e - →ϕ ϕ m Z ' = M ϕ χ ( Dashed ) m Z ' = M ϕ χ ( Dashed ) M ϕ , χ ( GeV ) C r o ss S e c t i on ( pb ) m Z ' = M ϕ χ ( Solid ) m Z ' = M ϕ χ ( Solid ) Δ ϕ χ = M ϕ χ , α D = ϵ = Δ ϕ χ = M ϕ χ , α D = ϵ = e + e - →χ χ γ e + e - →χ χ γ e + e - →ϕ ϕ γ e + e - →ϕ ϕ γ m Z ' = M ϕ χ ( Dashed ) m Z ' = M ϕ χ ( Dashed ) M ϕ , χ ( GeV ) C r o ss S e c t i on ( pb ) FIG. 1: The relations between M φ ,χ (GeV) and σ (pb) for e + e − → φ φ and e + e − → χ χ processes (left panel) and e + e − → φ φ γ and e + e − → χ χ γ processes (right panel). B. The numerical results
We first generate both scalar and fermion inelastic DM UFO model files from
Feyn-Rules [69], and then calculate cross sections of e + e − → φ φ ( χ χ ) and e + e − → φ φ ( χ χ ) γ processes via Madgraph5 aMC@NLO . The beam energies are set to be E ( e + ) = 4 . E ( e − ) = 7 . M φ ,χ and m Z (cid:48) , we study m Z (cid:48) = 2 . M φ ,χ and m Z (cid:48) = 3 M φ ,χ with fixed α D ≡ g D / π = 0 . (cid:15) = 0 .
01 and ∆ φ,χ = 0 . M φ ,χ . The relationbetween M φ ,χ (GeV) and σ (pb) for e + e − → φ φ ( χ χ ) processes are shown in the leftpanel of Fig. 1. Since all of these cross sections are proportional to (cid:15) α D , it is straightforwardto rescale cross sections with different values of (cid:15) and α D . On the other hand, the influencefrom the changes of ∆ φ,χ to cross sections are mild for ∆ φ,χ < . M φ ,χ . The scalar andfermion pair production cross sections can be scaled by β / and β / respectively, where β is the velocity of the final state particle in the CM frame. Because of this extra β factor forthe scalar case, cross sections for e + e − → φ φ are suppressed compared with e + e − → χ χ .On the other hand, the relation between M φ ,χ (GeV) and σ (pb) for e + e − → φ φ ( χ χ ) γ processes are shown in the right panel of Fig. 1 with the same parameter settings. Here, thebasic cuts E ( γ ) > . η < .
203 are applied for the ISR photon. We find the maincontributions of e + e − → φ φ ( χ χ ) γ come from e + e − → Z (cid:48) γ → φ φ ( χ χ ) γ processes. Itexplains why scalar and fermion pair cross sections in this process are very close to each The Z (cid:48) boson total decay width is automatically calculated in Madgraph5 aMC@NLO. /σ )( dσ/d cos θ ) distributions for e + e − → φ φ in CM frame (top left panel)and Belle II LAB frame (bottom left panel) and e + e − → χ χ in CM frame (top rightpanel) and Belle II LAB frame (bottom right panel).other in the right panel of Fig. 1.We then turn to the study of (1 /σ )( dσ/d cos θ ) distributions for e + e − → φ φ and e + e − → χ χ processes with fixed α D = 0 . (cid:15) = 0 .
01, ∆ φ,χ = 0 . M φ ,χ , and m Z (cid:48) = 3 M φ ,χ .The results in the CM and Belle II LAB frame are shown in Fig. 2. It is clear that the(1 /σ )( dσ/d cos θ ) distributions for e + e − → φ φ in the CM frame is proportional to 1 − cos θ and e + e − → χ χ behaves close to 1 + cos θ with some distortions from M χ , effects.However, since the measurement of time-of-flight for the long-lived particle is poor at theBelle II machine, we will lose this information and cannot make the Lorentz transformationfor the four-vector of displaced vertex from the Belle II LAB frame to the CM frame. Hence,the distributions in the CM frame is only a reference for the comparison with the ones inthe Belle II LAB frame. Hopefully, the real situation for the distributions in the Belle IILAB frame is just shifted to the electron beam direction for the Belle II machine comparedwith the ones in the CM frame. It is still clear to see the different distributions between12IG. 3: The (1 /σ )( dσ/d cos θ ) distributions for e + e − → φ φ γ in CM frame (top leftpanel) and Belle II LAB frame (bottom left panel) and e + e − → χ χ γ in CM frame (topright panel) and Belle II LAB frame (bottom right panel).scalar and fermion inelastic DM models.Finally, we fix the same parameter settings as Fig. 2 for the (1 /σ )( dσ/d cos θ ) distribu-tions of e + e − → φ φ ( χ χ ) γ processes in the CM and LAB frame in Fig. 3. Again, thedistributions in the CM frame is only for the comparison and it is clear to see the skew-ing behavior for the differential angular distributions of cross sections in the LAB framecompared with the ones in the CM frame. As discussed in the previous discussion, be-cause these distributions are highly dependent on the z parameter of the ISR photon, the e + e − → φ φ ( χ χ ) γ processes are not ideal for the DM spin discrimination. Only whenthe Z (cid:48) is on-shell produced and heavy, the ISR photon becomes soft ( z is small) such thatdifferences of (1 /σ )( dσ/d cos θ ) distributions from scalar and fermion inelastic DM modelsshow up.In summary, the size for cross sections of e + e − → φ φ ( χ χ ) and their (1 /σ )( dσ/d cos θ )distributions in the Belle II LAB frame can help us to statistically distinguish fermion13nelastic DM model from the scalar one even for the same parameter settings. IV. SEARCH FOR LONG-LIVED PARTICLES IN INELASTIC DARK MATTERMODELS AT BELLE II
In this section, we study how to search for long-lived particles in inelastic DM models atBelle II experiment. We first briefly overview the Belle II experiment and signal signatures,and then show some interesting kinematic distributions based on four signal benchmarkpoints. We further set up event selections for dilepton displaced signature and relevantresults are discussed, especially for the discrimination of scalar and fermion inelastic DMmodels. Finally, we solve kinematic equations of e + e − → χ χ → χ χ l + l − and e + e − → Z (cid:48) γ → χ χ γ → χ χ l + l − γ processes to determine both M χ and M χ . A. The Belle II experiment and signal signatures
The SuperKEKB accelerator of Belle II experiment is a circular asymmetric e + e − colliderwith the nominal collision energy of √ s = 10 .
58 GeV. The beam parameters are E ( e + ) = 4GeV and E ( e − ) = 7 GeV. The planned full integrated luminosity for the final dataset is50 ab − . In this work, the following Belle II sub-detectors are relevant : the tracking systemincluding vertex detectors (VXD) and central drift chamber (CDC), the electromagneticcalorimeter (ECAL), and muon system.We consider the following two kinds of processes in inelastic DM models : e + e − → φ φ ( χ χ ) → φ φ ( χ χ ) l + l − , (28) e + e − → Z (cid:48) γ → φ φ ( χ χ ) γ → φ φ ( χ χ ) l + l − γ . (29)And there are five possible signatures can be searched for at the Belle II experiment : • Mono- γ : If the excited DM state decays outside the detector or the decay productsare too soft to be detected in Eq.(29). • Mono- γ with prompt lepton pair: The excited DM state promptly decays and thevisible products can be successfully detected and defined in Eq.(29). Because the major contribuitons of e + e − → φ φ ( χ χ ) γ come from the on-shell Z (cid:48) production, hence, weonly focus on the process e + e − → Z (cid:48) γ → φ φ ( χ χ ) γ in the following analysis unless noted otherwise. Mono- γ with displaced lepton pair : The excited DM state is long-lived and decaysinside the detector leaving the displaced vertex in Eq.(29). • Only prompt lepton pair : The same as the second one but without ISR photon inEq.(28). • Only displaced lepton pair : The same as the third one but without ISR photon inEq.(28).The Mono- γ signature is well-studied as shown in Ref. [21–23]. Since we focus on searchingfor long-lived particles in this paper, only signatures with displaced vertex are studied. Onthe other hand, the signature with displaced charged pion pair is also possible from thelong-lived excited DM state decay. However, for simplicity, here we only study signatureswith displaced electron and muon pair.In order to make our results more realistic, we follow Ref. [11] to involve the detectorresolution effects with Gaussian smearing at Belle II experiment. The tracking resolution ofelectron/muon momentum in the CDC is given by σ p lt /p lt = 0 . p lt [GeV] ⊕ . /β, (30)where p lt is the transverse momentum of electron/muon track and β is its velocity in thenatural unit. We conservatively apply σ p l /p l = 0 .
005 in our event analysis. On the otherhand, the muon efficiency in the muon system is approximated to 0 .
98. The photon mo-mentum resolution in the ECAL is approximated to σ E γ /E γ = 2% where E γ is the energyof photon. Finally, we use the resolution of σ r DV = 26 µ m for the displaced vertex vector oflepton pair in our analysis. B. Benchmark points and kinematic distributions
In this study, we use the following four benchmark points (BPs) to display our analysis : • (I) M φ ,χ = 0 . φ ,χ = 0 . M φ ,χ , m Z (cid:48) = 3 M φ ,χ and (cid:15) = 2 × − • (II) M φ ,χ = 3 . φ ,χ = 0 . M φ ,χ , m Z (cid:48) = 3 M φ ,χ and (cid:15) = 2 × − • (III) M φ ,χ = 1 . φ ,χ = 0 . M φ ,χ , m Z (cid:48) = 2 . M φ ,χ and (cid:15) = 10 − e + e − → χ χ → χ χ l + l − based on the fourBPs in the main text. Top left panel: lepton pair invariant mass, M ll (GeV), top rightpanel: projection of the decay vertex distance on z axis, r z (mm), bottom left panel: leptonpair energy, E ll (GeV), bottom right panel: angular distance between lepton pair, ∆ R ll . • (IV) M φ ,χ = 2 . φ ,χ = 0 . M φ ,χ , m Z (cid:48) = 2 . M φ ,χ and (cid:15) = 10 − with fixed α D = 0 .
1. The first BP is inspired from the allowed parameter space to explainthe muon anomalous magnetic moment excess [24]. Other three BPs are the ones whichcannot be covered from Mono- γ searches from BaBar and Belle II [21]. In addition, theycan also be searched for in the future long-lived particle experiments, like CODEX-b [70],SeaQuest [71], FASER [72], and MATHUSLA [73].We show some interesting kinematic distributions based on the above four BPs for e + e − → χ χ → χ χ l + l − and e + e − → Z (cid:48) γ → χ χ γ → χ χ l + l − γ in Fig. 4 and Fig. 5, respectively.First, the lepton pair invariant mass is smaller or equal to the mass splitting ∆ χ (dashedlines), 2 m l ≤ M ll ≤ ∆ χ . Once we have enough signal events, the threshold value of M ll canhelp us to roughly determine the mass splitting ∆ χ in inelastic DM models. Second, thedistance of displaced vertex of χ is not only dependent on M χ , ∆ χ , (cid:15) and α D , but also theboost of χ in the LAB frame. Here we show projection of the decay vertex distance on z axis, r z . It is clear to see that the χ decay length in BP1 (BP2) is the shortest (longest) one16IG. 5: Various kinematic distributions for e + e − → Z (cid:48) γ → χ χ γ → χ χ l + l − γ based onthe four BPs in the main text. Top left panel: lepton pair energy, E ll (GeV), top rightpanel: angular distance between lepton pair, ∆ R ll , middle left panel: projection of thedecay vertex distance on z axis, r z (mm), middle right panel: the maximum angulardistance between ISR photon and one of the lepton in lepton pair, ∆ R maxγl , bottom leftpanel: photon energy in the CM frame, E ( γ ) (GeV), bottom right panel : reconstructed Z (cid:48) invariant mass, m Z (cid:48) (GeV).on average. The distributions for projection of the decay vertex distance on the transverseplane, R xy , are similar. On the other hand, E ll and ∆ R ll are energy and angular distance forthe lepton pair. We can find the boost of χ is smaller for the process involving ISR photonby comparing r z , E ll and ∆ R ll distributions in Fig. 4 and Fig. 5. The larger boost on χ causes the longer r z , larger E ll , and smaller ∆ R ll distributions. Besides, we will see larger E ll ,and smaller r z distributions of BP1 suffer from more severe smearing effects from detector17esolutions such that it is more challenge to reconstruct DM mass and mass splitting betweenDM excited and ground states in the real experiment. Third, if the Z (cid:48) is light enough, theISR photon and Z (cid:48) are almost back-to-back produced such that the lepton from Z (cid:48) decaycan have large angular distance with the ISR photon. Here we show the maximum angulardistance between ISR photon and one of the lepton in lepton pair, ∆ R maxγl . Fourth, we boostthe E ( γ ) distribution of ISR photon from the Belle II LAB frame to the CM frame and itis clear to see its mono-energetic behavior. This feature can largely reduce SM backgroundswhich are continuous in the photon energy spectrum. Finally, the Z (cid:48) invariant mass can bereconstructed from Eq.(39) with measurable observables. Because of the detector resolution,it is clear to find smearing effects from the true Z (cid:48) mass (dashed lines), especially for the caseof light Z (cid:48) . We will see this phenomenon critically affects abilities to determine DM massand the mass splitting between DM excited and ground states for the process in Eq.(29).Notice those kinematic distributions for the scalar inelastic DM model are similar to Fig. 4and Fig. 5, so we do not show them again. C. Event selections and results
We closely follow Ref. [21] to set up event selections for the displaced signature at BelleII. We only conservatively consider the following two background-free regions after eventselections in our analysis : • low R xy region (100% detection efficiency) : 0 . < R xy ≤ . . 03 GeV,and their opening angle larger than 0 . . < R xy ≤ . R xy region can be extended.We summarize all event selections in our analysis in Table I. For the muon pair in the finalstate, we veto the invariant mass region, 0 . ≤ m µ + µ − ≤ . 52 GeV, to reject backgroundsfrom K S decay. More discussions on the trigger issue can be found in Ref. [21]. Here, event18ABLE I: Event selections for the displaced vertex analysis at Belle II. Objects Selectionsdisplaced vertex (i) − 55 cm ≤ z ≤ 140 cm(ii) 17 ◦ ≤ θ DVLAB ≤ ◦ electrons (i) both E ( e + ) and E ( e − ) > . θ ee > . m ee > . 03 GeVmuons (i) both p T ( µ + ) and p T ( µ − ) > . 05 GeV(ii) opening angle of pair θ µµ > . m µµ > . 03 GeV(iv) veto 0 . 48 GeV ≤ m µµ ≤ . 52 GeVphotons (i) E γ LAB > . ◦ ≤ θ γ LAB ≤ ◦ selections in Table I are done offline and we simply assume these events are already triggeredand stored.The results for four BPs are shown in Table II and III for processes in Eq.(28) and (29),respectively. Eff.(low R xy ) and Eff.(high R xy ) are efficiencies of low R xy and high R xy regions after involving the event selections. Here we use an integrated luminosity of 1 ab − to calculate number of signal events ( N event ). We can find most events in BP1 and BP3(BP2 and BP4) are located in low (high) R xy regions.We further compare polar angular distributions in the Belle II LAB frame after eventselections in Table I for processes in Eq.(28) and (29) with electron pair in the final statein Fig. 6 and 7 with N event in Table II and III, separately. In Fig. 6, since the selectedregions can be approximately background free, we can directly count the events bin by binto distinguish the signals are either from scalar or fermion inelastic DM models for theabove four BPs even only L = 1 ab − is used. However, in Fig. 7, it is hard to distinguishthe signals are either from scalar or fermion inelastic DM models except for the BP2. Wehave discussed this behavior in Sec. III that we only have the chance to separate these twomodels if the ISR photon is soft enough. Therefore, in addition to the usual search processin Eq.(29) at B-factories, we propose to consider the process in Eq.(28) with the displaced19IG. 6: The polar angular distributions in the Belle II LAB frame after event selections inTable I for process in Eq.(28) with electron pair in the final state and N event in Table II.FIG. 7: The same as Fig. 6, but for process in Eq.(29) with electron pair in the final stateand N event in Table III.20ABLE II: The analysis results of four BPs for process in Eq.(28). The first column is thetype of inelastic DM models. The second column is labeled by BP which e and µ denotethe type of lepton pair in the final state. The third column is production cross sections forprocess in Eq.(28). The fourth and fifth columns are efficiencies of low R xy and high R xy regions defined in the main text, respectively. Final column is the number of signal eventswith an integrated luminosity of 1 ab − . Type BP σ (fb) Eff.(low R xy ) Eff.(high R xy ) N event scalar BP1 e . 14 16 . 98 0% 1 . × BP2 e . 39 0 . 15% 2 . 48% 1 . × BP2 µ . 15 0 . 21% 3 . 33% 217 . e . 86 10 . 06% 0 . 70% 200 . µ . 61 11 . 25% 0 . 74% 73 . e . 23 1 . 56% 9 . 34% 243 . µ . 74 1 . 72% 10 . 78% 92 . e . 00 14 . 26% 0% 5 . × BP2 e . 80 0 . 17% 2 . 35% 1 . × BP2 µ . 63 0 . 22% 2 . 97% 1 . × BP3 e . 99 10 . 20% 0 . 42% 848 . µ . 69 11 . 20% 0 . 46% 313 . e . 71 1 . 57% 7 . 82% 1 . × BP4 µ . 88 1 . 69% 8 . 75% 405 . vertex trigger suggested in Ref. [21] which is more sensitive to distinguish the spin of DMin inelastic DM models.Finally, we use the most aggressive integrated luminosity with 50 ab − for event selectionsin Table I at Belle II for scalar and fermion inelastic DM models and predict the futurebounds from e + e − → φ φ ( χ χ ) and e + e − → φ φ ( χ χ ) γ processes in Fig. 8. Notice theoff-shell Z (cid:48) production in the second process is also considered. Here we fix the parameters, α D = 0 . m Z (cid:48) = 3 M φ ,χ and ∆ φ,χ = 0 . M φ ,χ , and apply 90% C.L. contours whichcorrespond to an upper limit of 2.3 events with the assumption of background-free. Notice21ABLE III: The same as Table II, but for process in Eq.(29). Type BP σ (fb) Eff.(low R xy ) Eff.(high R xy ) N event scalar BP1 e . 70 6 . 70% 0% 1 . × BP2 e . 85 0 . 16% 2 . 27% 3 . × BP2 µ . 85 0 . 20% 2 . 87% 517 . e . 13 7 . 64% 0 . 02% 392 . µ . 69 8 . 83% 0 . 03% 149 . e . 14 1 . 86% 3 . 29% 367 . µ . 35 2 . 02% 2 . 87% 114 . e . 60 6 . 14% 0% 1 . × BP2 e . 10 0 . 16% 2 . 16% 3 . × BP2 µ . 66 0 . 18% 2 . 67% 503 . e . 05 7 . 77% 0 . 02% 393 . µ . 70 8 . 89% 0 . 02% 151 . e . 14 1 . 95% 3 . 14% 363 . µ . 37 2 . 05% 3 . 44% 130 . constraints from model-independent LEP bound [74] and BaBar mono- γ bound [23] areadded for the comparison. We closely follow Ref. [21] for the recasting of BaBar mono- γ constraints in inelastic DM models. For M φ ,χ < . M φ ,χ (cid:38) . g − µ [24], we perform the same22 aBar mono - γ BaBar mono - γ e + e - →ϕ ϕ γ e + e - →ϕ ϕ γ e + e - →ϕ ϕ e + e - →ϕ ϕ LEPLEP - - M ϕ ( GeV ) ϵ α D = m Z ' = M ϕ , Δ ϕ = M ϕ BaBar mono - γ BaBar mono - γ e + e - →χ χ γ e + e - →χ χ γ e + e - →χ χ e + e - →χ χ LEPLEP - - M χ ( GeV ) ϵ α D = m Z ' = M χ , Δ χ = M χ FIG. 8: The future bounds from e + e − → φ φ ( χ χ ) and e + e − → φ φ ( χ χ ) γ processesfor event selections in Table I with the integrated luminosity of 50 ab − . Here parameters α D = 0 . m Z (cid:48) = 3 M φ ,χ and ∆ φ,χ = 0 . M φ ,χ are fixed and 90% C.L. contours whichcorrespond to an upper limit of 2.3 events with the assumption of background-free areapplied. The model-independent LEP bound [74] and BaBar mono- γ bound [23] are alsoshown. - - M ϕ ( GeV ) ϵ α D = m Z ' = M ϕ , Δ ϕ = M ϕ LEP BaBar mono - γ σ f a v o r e d a μ σ e x c l u d e d a μ e + e - →ϕ ϕ γ e + e - →ϕ ϕ - - M χ ( GeV ) ϵ α D = m Z ' = M χ , Δ χ = M χ LEP BaBar mono - γ σ f a v o r e d a μ σ e x c l u d e d a μ e + e - →χ χ e + e - →χ χ γ FIG. 9: The same as Fig. 8 but for m Z (cid:48) = 3 M φ ,χ and ∆ φ,χ = 0 . M φ ,χ . The greenshaded region bounded by the green dashed lines is the 2 σ allowed region for the ( g − µ excess and the lighter gray region excluded by the ( g − µ at 5 σ C.L.analysis in Fig. 9 for α D = 0 . m Z (cid:48) = 3 M φ ,χ and ∆ φ,χ = 0 . M φ ,χ fixed. We show the 2 σ allowed and the 5 σ excluded regions for the muon ( g − µ as well as the model-independentLEP bound [74] and the BaBar mono- γ bound [23]. Once the mass of DM is increasing,the mass splitting between DM ground and excited states is also enhanced, hence, the final23tate lepton pair becomes more energetic. Therefore, one can find the BaBar mono- γ boundis slightly weakened when the displaced channel is open. Since the recasting in Ref [24] doesnot include the requirement on the angle θ DVLAB in Table I, the results are weaker than theones in Ref. [21]. That is the reason why this parameter space is still allowed in Ref [24].In order to make a more precise recasting, we follow Ref. [21] for the BaBar mono- γ boundand it is shown in dotted gray line in Fig. 9 which can already cover the ( g − µ allowedwindow in the parameter space. Even the recasting BaBar mono- γ bound can close thisregion, we find our result can give much stronger bound, (cid:15) ∼ O (10 − − − ), especially inthe parameter space where the mono- γ searches are suppressed. Therefore we can explicitlycover this area using displaced lepton search at the Belle II in a very first stage of the run.We found the future bounds from the process without ISR photon can be stronger than theone with ISR photon in Fig. 8. On the other hand, in the lower M φ region, the boundfor process with ISR photon get stronger because its cross sections are larger than the oneswithout ISR photon as shown in Fig. 1. D. Determination of the DM mass and mass splitting of DM sector In general, we cannot uniquely determine the DM mass at colliders because of lackingenough constraints for invisible particles in the final state. We can take a mono- γ searchfor DM pair production at B-factories as an example. If the DM is denoted by χ , theprocess for mono- γ search at B-factories is e + e − → χχγ (or e + e − → χχγ ). There are eightunknown values from four-momentum of two DMs in the final state. Unfortunately, onlyfive constraints in this process : four from the four-momentum conservation and one fromthe same mass for the DM pair. We still need three extra conditions to uniquely determinethe DM mass for each event.Now we turn our attention to the case of inelastic DM models. For the process in Eq.(28)or (29), if the φ ( χ ) is long-lived and leave the displaced vertex at the Belle II detectors,we will have two extra constraints. Again, there are still eight unknown values from four-momentum of two φ (cid:48) s ( χ (cid:48) s ) in the final state. However, because of the charge neutrality ofthe φ ( χ ), a three-momentum vector of φ ( χ ) is proportional to the direction of displacedvertex (DV) [57, 60] −→ p φ ( χ ) = |−→ p φ ( χ ) | (cid:98) r DV , (31)24here −→ p φ ( χ ) is the three-momentum vector of φ ( χ ) and (cid:98) r DV is the unit vector of displacedvertex from φ ( χ ). Therefore, we have two more constraints, and there are seven constraintsfor this kind of processes in total. We still need one more condition to uniquely determinethe DM mass. Let’s first study the kinematic equation for processes in Eq.(28) and (29)and solve them event-by-event from our Monte Carlo samples and then involve the detectorresolution effects .For e + e − → χ χ → χ χ l + l − , we first set the direction of DV and following fourmomentum of e + , e − , χ and χ as (cid:98) r DV = (sin θ cos φ, sin θ sin φ, cos θ ) , (32)and p e + = ( E + , , , E + ) ,p e − = ( E − , , , − E − ) ,p χ = ( E χ , |−→ p χ | (cid:98) r DV ) ,p χ = ( E χ , −→ p χ ) , (33)with E − (cid:54) = E + in the Belle II LAB frame, then according to energy and momentum conser-vation, we can write down E χ = E + + E − − E χ , −→ p χ = ( −|−→ p χ | sin θ cos φ, −|−→ p χ | sin θ sin φ, E + − E − − |−→ p χ | cos θ ) . (34)After applying the above two equations to E χ i = |−→ p χ i | + M χ i for i = 1 , 2, the E χ can bewritten as E χ = 12 (cid:2) sin θ ( E + E − ) + 2(1 + cos θ ) E + E − (cid:3) [( E + + E − )(4 E + E − + M χ − M χ ) ± | ( E − − E + ) cos θ | (( M χ − E + E − ) − θ ( E + E − ) + 4 cos θE + E − + M χ ] M χ + M χ ) / ] . (35)We then consider the energy and momentum conservation for χ → χ f f , E χ = E χ (cid:48) + E f + E f ≡ E χ (cid:48) + E V (cid:48) , |−→ p χ | (cid:98) r DV = −→ p χ (cid:48) + −→ p f + −→ p f ≡ −→ p χ (cid:48) + −→ p V (cid:48) . (36) Here we use the fermion inelastic DM model as an example. The same method can be applied to thescalar inelastic DM model, so we will not repeat it thereafter. .200 0.225 0.250 0.275 0.300 0.325 0.350 0.375 0.400 M (GeV)0.100.110.120.130.140.15 ( G e V ) M maxl + l = event-1event-2event-3BP1 : e + e l + l M (GeV)0.2000.2250.2500.2750.3000.3250.3500.3750.400 ( G e V ) M maxl + l = event-1 event-2event-3 BP2 : e + e l + l M (GeV)0.3000.3250.3500.3750.4000.4250.4500.4750.500 ( G e V ) M maxl + l = event-1event-2 event-3BP3 : e + e l + l M (GeV)0.3000.3250.3500.3750.4000.4250.4500.4750.500 ( G e V ) M maxl + l = event-1event-2event-3 BP4 : e + e l + l FIG. 10: The solutions for Eq.(37) with E χ in Eq.(35) of e + e − → χ χ → χ χ l + l − process on the ( M χ , ∆ χ ) plane for four BPs. Here we display three arbitrary Monte Carloevents for each BP with red, green and blue lines. On the other hand, the black line showsthe kinematic endpoint measurement of M maxl + l − .Finally, we can receive the following equation for inputs of E V (cid:48) , −→ p V (cid:48) and (cid:98) r DV , M χ − M χ − E χ E V (cid:48) + E V (cid:48) − |−→ p V (cid:48) | + 2 (cid:113) E χ − M χ ( (cid:98) r DV · −→ p V (cid:48) ) = 0 , (37)where E χ is shown in Eq.(35). In Fig. 10, we display three arbitrary Monte Carlo eventson the ( M χ , ∆ χ ) plane for four BPs according to Eq.(37). On the other hand, the blackline shows the kinematic endpoint measurement M maxl + l − . We can find that the lines from allevents and M maxl + l − cross to the same point which is the true ( M χ , ∆ χ ) in our four BPs. Thisis a simple application of the Kinematic Focus Point Method proposed in Ref. [58].Similarly, for e + e − → Z (cid:48) γ → χ χ γ → χ χ l + l − γ , we use the same settings inEq.(32),(33) and set the four momentum of ISR photon in the Belle II LAB frame as p γ = ( E γ , −→ p γ ) . (38)According to energy and momentum conservation, the four momentum of on-shell Z (cid:48) can26 .200 0.225 0.250 0.275 0.300 0.325 0.350 0.375 0.400 M (GeV)0.100.110.120.130.140.15 ( G e V ) M maxl + l = event-1 event-2event-3BP1 : e + e Z l + l M (GeV)0.2000.2250.2500.2750.3000.3250.3500.3750.400 ( G e V ) M maxl + l = event-2event-3event-1BP2 : e + e Z l + l M (GeV)0.300.350.400.450.500.550.600.650.70 ( G e V ) M maxl + l = event-1event-2event-3BP3 : e + e Z l + l M (GeV)0.3000.3250.3500.3750.4000.4250.4500.4750.500 ( G e V ) M maxl + l = event-3event-2event-1BP4 : e + e Z l + l FIG. 11: The same as Fig. 10, but the solutions for Eq.(37) with E χ in Eq.(41) of e + e − → Z (cid:48) γ → χ χ γ → χ χ l + l − γ process.be written as p Z (cid:48) = ( E Z (cid:48) , −→ p Z (cid:48) ) = ( E + + E − − E γ , − p x,γ , − p y,γ , E + − E − − p z,γ ) , (39)where −→ p γ = ( p x,γ , p y,γ , p z,γ ) and p Z (cid:48) = m Z (cid:48) . Furthermore, we represent the four momentumof χ as p χ = ( E Z (cid:48) − E χ , −→ p Z (cid:48) − |−→ p χ | (cid:98) r DV ) , (40)and write down the E χ as E χ = 12 [ E Z (cid:48) − ( (cid:98) r DV · −→ p Z (cid:48) ) ] [ E Z (cid:48) ( m Z (cid:48) + M χ − M χ ) ± | (cid:98) r DV · −→ p Z (cid:48) | (cid:113) ( m Z (cid:48) + M χ − M χ ) − M χ [ E Z (cid:48) − ( (cid:98) r DV · −→ p Z (cid:48) ) ]] . (41)Finally, the kinematic equation for χ → χ f f is the same as Eq.(37) with E χ in Eq.(41).Again, we show three arbitrary Monte Carlo events on the ( M χ , ∆ χ ) plane for four BPsaccording to Eq.(37) in Fig. 11. All events and M maxl + l − cross to the true ( M χ , ∆ χ ) in ourfour BPs. 27IG. 12: The solutions for Eq.(37) with E χ in Eq.(35) of e + e − → χ χ → χ χ l + l − process on the ( M χ , M χ ) plane for four BPs after involving detector resolution effects andevent selections. Here 100 signal events are taken into account for each BP and the bin sizefor both x and y axes are set to be 10 MeV.In order to make the Kinematic Focus Point Method fit to reality, we involve the detectorresolution effects from Sec. IV A and event selections in Table. I. Here we conservativelyassume there are 100 signal events for each BP can be recorded. Since we can solve theabove kinematic equations for any two signal events, we will get C = 4950 solutions from100 signal events. After removing the unphysical solutions, we show the distributions ofthese solutions on the ( M χ , M χ ) plane for four BPs of e + e − → χ χ → χ χ l + l − and e + e − → Z (cid:48) γ → χ χ γ → χ χ l + l − γ processes in Figs. 12 and 13, respectively. Here the binsize for both x and y axes are set to be 10 MeV. Moreover, the most probable mass valuesand their statistical errors for four BPs in our simulation are summarized in Table. IV and V.The statistical errors are estimated as the root mean square (rms) value with respective tothe most probable values rms = (cid:113)(cid:80) N phys i =1 ( M i − M peak ) /N phys , where N phys is the number28IG. 13: The same as Fig. 12, but the solutions for Eq.(37) with E χ in Eq.(41) of e + e − → Z (cid:48) γ → χ χ γ → χ χ l + l − γ process.of physical solutions from 100 signal events and M peak is the most probable value.In Fig. 12, except for BP1 on top left panel, M χ and M χ can be determined within thepercentage of deviation for other three BPs. As shown in Fig. 4, both the decay length of χ is relative short and E ll is relative large for BP1 compared with other three BPs such thatevents in BP1 cannot be avoided to severely suffer from detector resolution effects. Thisexplains why the physical solutions are reduced and the reconstructed M χ and M χ areshifted to larger values with higher smearing for BP1. In Fig. 13, M χ and M χ can stillbe determined within the percentage of deviation for BP2 and BP4. However, apart fromdetector resolution effects from charged leptons and displace vertex, the ability to preciselypin down M χ and M χ for the process in Eq.(29) also relies on how well the on-shell Z (cid:48) can be reconstructed. As shown in Fig. 5, the peak of m Z (cid:48) is rather spread out (or evenshifted) for BP1 and BP3 than for BP2 and BP4. According to Eq.(39), the reconstructed29 P N phys ( M χ , M χ ) true rms( M χ , M χ ) peak BP1 4473 (0 . , . 30) (0 . , . . , . . , . 00) (0 . , . . , . . , . 00) (0 . , . . , . . , . 00) (0 . , . . , . TABLE IV: The peak measured values and root mean square (rms) for ( M χ , M χ ) of fourbenchmark points in e + e − → χ χ → χ χ l + l − process. N phys is the number of physicalsolutions from 100 signal events. All numbers are in GeV unit. BP N phys ( M χ , M χ ) true rms( M χ , M χ ) peak BP1 901 (0 . , . 30) (0 . , . . , . . , . 00) (0 . , . . , . . , . 00) (0 . , . . , . . , . 00) (0 . , . . , . TABLE V: The same as Table. IV, but in e + e − → Z (cid:48) γ → χ χ γ → χ χ l + l − γ process. Z (cid:48) four-momentum is highly dependent on the detector resolution effects of ISR photon.Once the E ( γ ) increases, its smearing is also enhanced such that the reconstruction of theon-shell Z (cid:48) becomes poor. In the end, the physical solutions of BP1 and BP3 are largely30educed and shifting behaviors of M χ and M χ become severe in Fig. 13 compared with theones in Fig. 12. This is another reason to encourage our experimental colleagues to searchfor not only the usual process in Eq.(29) at B-factories, but also the process in Eq.(28) withthe displaced vertex trigger which is more sensitive to determine the mass and mass splittingof DM sector for some parameter space in the inelastic DM models. V. CONCLUSION The dark sectors with light DM candidate are still all the rage as the null signal resultfrom DM direct detetion has been reported until now. On the other hand, the quietnessof beyond the SM searches at the LHC also hints new physics signatures may hide in theelusive corner. New kind of signatures such as long-lived particles at colliders become moreand more popular and may guide a royal road for new physics evidences. Among thesedark sector models, the inelastic DM model is an appealing example which can both allowlight DM candidate and predict long-lived particle which can be searched for at colliders.Therefore, we focus on the displaced vertex signatures in inelastic DM models at Belle II inthis work.In the inelastic DM models, if the mass splitting between DM excited and ground statesis small enough, the co-annihilation becomes the dominant channel for thermal relic densityand the DM excited state can be long-lived at the collider scale. We first review scalar andfermion inelastic DM models with U (1) D gauge symmetry in Sec. II and point out the darkHiggs sector caused the mass splitting between DM excited and ground states and provideda mechanism to the Z (cid:48) mass. Besides, the off-diagonal interaction between Z (cid:48) and DM sectoris derived from the covariant derivative term of DM sector.The analytical representations and numerical results for cross sections of e + e − → φ φ ( χ χ ) and e + e − → φ φ ( χ χ ) γ processes are studied in Sec. III. In the first pro-cess, there is an extra β factor suppression for boson pair cross sections compared withfermion pair ones. In the second process, the dominant channel comes from e + e − → Z (cid:48) γ → φ φ ( χ χ ) γ . Therefore, cross sections in this process for scalar and fermion inelatic DMmodels are very close to each other. In addition, the polar angle distributions of φ ( χ ) atthe Belle II LAB frame as shown in Fig. 2 are helpful to distinguish these two models.The novel dilepton displaced vertex signatures at Belle II are the main targets of this31tudy. We include detector resolution effects in Sec. IV A and event selections in Table. I forsignal processes in Eq.(28) and (29). In Fig. 6, we display four BPs of scalar and fermioninelastic DM models and these two models can be discriminated even with an integratedluminosity less than 1 ab − . Our analysis results are summarized in Figs. 8 and 9, whichindicate that the future bounds from dilepton displaced vertex searches in inelastic DMmodels at Belle II are stronger than previous constraints for two benchmark mass windows0 . (cid:46) M φ ,χ (cid:46) . . (cid:46) M φ ,χ (cid:46) . g − µ . Thereforethe early stage of the Belle II experiment can explicitly close off this area in the parameterspace. We further apply the Kinematic Focus Point method for e + e − → χ χ → χ χ l + l − and e + e − → Z (cid:48) γ → χ χ γ → χ χ l + l − γ processes to determine both M χ and M χ masses.Especially, we find the mass and mass splitting of DM sector can be determined within thepercentage of deviation as shown in Table. IV and V. ACKNOWLEDGMENT We thank Myeonghun Park and Youngjoon Kwon for very helpful discussions and use-ful information. The work is supported in part by KIAS Individual Grants, Grant No.PG076201 (DWK), PG021403 (PK) and PG075301 (CTL) at Korea Institute for Ad-vanced Study, and by National Research Foundation of Korea (NRF) Grant No. NRF-2019R1A2C3005009 (PK), funded by the Korean government (MSIT).32 ppendix A: The analytical representations for cross sections of e + e − → φ φ ( χ χ ) in the LAB frame In this appendix, we show the analytical representations for cross sections of e + e − → φ φ ( χ χ ) in the LAB frame which have been mentioned in Sec. III. We first set the fourmomentum of the initial e + , e − and φ ( χ ), φ ( χ ) in the LAB frame as p e + = ( E + , , , E + ) ,p e − = ( E − , , , − E − ) ,p φ ( χ ) = ( E φ ( χ ) , p x,φ ( χ ) , p y,φ ( χ ) , p z,φ ( χ ) ) ,p φ ( χ ) = ( E φ ( χ ) , p x,φ ( χ ) , p y,φ ( χ ) , p z,φ ( χ ) ) . (A1)According to the four momentum conservation for this scattering process, p φ ( χ ) can bewritten as p φ ( χ ) = ( E + + E − − E φ ( χ ) , − p x,φ ( χ ) , − p y,φ ( χ ) , E + − E − − p z,φ ( χ ) ) . (A2)The dσ/d cos θ for e + e − → φ φ via s-channel Z (cid:48) in the LAB frame can be written as dσ ( e + e − → φ φ ) d cos θ (cid:12)(cid:12)(cid:12)(cid:12) LAB = (cid:15) e g D πE + E − [(4 E + E − − m Z (cid:48) ) + m Z (cid:48) Γ Z (cid:48) ] × [8 E + E − ( E φ ( E + + E − ) + ( E − − E + ) |−→ p φ | cos θ ) − E + E − ( M φ + 3 M φ ) + ( M φ − M φ ) + (cid:0) E + ( E φ − E − − |−→ p φ | cos θ ) + M φ − M φ (cid:1) × (cid:0) E − ( E φ − E + + |−→ p φ | cos θ ) + M φ − M φ (cid:1) ] × (cid:18) − dE φ d cos θ + cos θ d |−→ p φ | d cos θ − sin θ |−→ p φ | (cid:19) , (A3)where |−→ p φ | = (cid:113) E φ − M φ and E φ is the same as Eq.(35) with χ , → φ , .On the other hand, the dσ/d cos θ for e + e − → χ χ via s-channel Z (cid:48) in the LAB frame33an be written as dσ ( e + e − → χ χ ) d cos θ (cid:12)(cid:12)(cid:12)(cid:12) LAB = (cid:15) e g D πE + E − [(4 E + E − − m Z (cid:48) ) + m Z (cid:48) Γ Z (cid:48) ] × [2( E + E − ) E χ (1 + cos θ ) − ( E − − E + )( M χ − M χ ) |−→ p χ | cos θ + ( E + + E − ) (cid:0) E − − E + ) |−→ p χ | cos θ + M χ − M χ (cid:1) E χ − E + E − ) M χ cos θ + 4 E + E − M χ M χ ] × (cid:18) − dE φ d cos θ + cos θ d |−→ p φ | d cos θ − sin θ |−→ p φ | (cid:19) , (A4)where |−→ p χ | = (cid:112) E χ − M χ and E χ can be found in Eq. (35). Appendix B: The full analytical representations of e + e − → φ φ ( χ χ ) → φ φ ( χ χ ) f f In this appendix, we display the full analytical representations of e + e − → φ φ ( χ χ ) → φ φ ( χ χ ) f f as mentioned in Sec. III. The differential cross section for e + ( p ) e − ( p ) → φ ( p ) φ ( p ) f ( p ) f ( p ) can be represented as dσ ( e + e − → φ φ f f ) = (2 π ) |M → | s × d Φ ( p + p ; p , p , p , p ) , (B1)where M → is the scattering amplitude for e + e − → φ φ f f , s is the square of centre-of-mass energy and Φ ( p + p ; p , p , p , p ) is the four-body phase space. Since φ is on-shell produced in this process, we apply the narrow width approximation for on-shell φ toseparate e + e − → φ φ f f to e + e − → φ φ with φ → φ f f . The differential phase space d Φ ( p + p ; p , p , p , p ) can be expanded as d Φ ( p + p ; p , p , p , p ) = d Φ ( p + p ; p , l ) × d Φ ( l ; p , p , p )(2 π ) dl , (B2)where l is the four momentum of on-shell φ . Therefore, we can disassemble the differentialcross section in the following form, dσ ( e + e − → φ φ f f ) = dσ ( e + e − → φ φ ) dB ( φ → φ f f ) , (B3)where dσ ( e + e − → φ φ ) = (2 π ) |M → | s × d Φ ( p + p ; p , l ) , (B4)34nd dB ( φ → φ f f ) = 2 π Γ φ δ ( l − M φ ) 12 M φ |M → | d Φ ( l ; p , p , p )(2 π ) dl = d Γ φ → φ ff Γ φ δ ( l − M φ )(2 π ) dl , (B5)where Γ φ is the total width of φ , M → is the scattering amplitude for e + e − → φ φ and M → is the decay amplitude for φ → φ f f . After applying the narrow width approxima-tion for on-shell φ , |M → | can be transferred to |M → | = |M → | × πδ ( l − M φ ) M φ Γ φ |M → | . (B6)Furthermore, we integrate out some phase space variables and leave the following kine-matic variables : • cos θ : θ is the direction of φ relative to the positron beam direction. • x − : x − = 2 E f /M φ in the φ rest frame. • x + : x + = 2 E f /M φ in the φ rest frame.for the representation of differential cross section, dσ ( e + e − → φ φ f f ) d cos θdx − dx + = dσ ( e + e − → φ φ ) d cos θ dB ( φ → φ f f ) dx − dx + . (B7)Notice the opening angle between the fermion pair can be written ascos θ ff = 1 − x − + x + − M φ /M φ ) / ( x − x + ) . (B8)The dσ ( e + e − → φ φ ) /d cos θ can be found in Eq.(22) in the CM frame and Eq.(A3) in theLAB frame and d Γ( φ → φ f f ) /dx − /dx + can be written as d Γ( φ → φ f f ) dx − dx + = g D (cid:15) αM φ π (1 − x + )(1 − x − ) M φ − ( M φ − ∆ φ ) (cid:2) M φ ( x + + x − − − m Z (cid:48) (cid:3) + m Z (cid:48) Γ Z (cid:48) . (B9)One can find x + and x − are symmetric in this formula.The associated formulae for e + ( p ) e − ( p ) → χ ( p ) χ ( p ) f ( p ) f ( p ) can be obtained fromEq.(B1) to (B7) by replacing φ , to χ , . The dσ ( e + e − → χ χ ) /d cos θ can be found in35q.(25) in the CM frame and Eq.(A4) in the LAB frame and d Γ( χ → χ f f ) /dx − /dx + canbe written as d Γ( χ → χ f f ) dx − dx + = g D (cid:15) αM χ π (cid:2) M χ ( x + + x − − − m Z (cid:48) (cid:3) + m Z (cid:48) Γ Z (cid:48) × {− M χ (cid:2) x + x − + 2( x + + x − ) (cid:3) + 4 M χ ∆ χ (1 + x + + x − ) − M χ (cid:2) (6 + x + + x − )∆ χ + 2(6 + x + + x − ) m f (cid:3) + 2 M χ (∆ χ + 4∆ χ m f ) − m f (∆ χ + 4 m f ) } . (B10) [1] C. S. Frenk and S. D. M. White, Annalen Phys. , 507-534 (2012)doi:10.1002/andp.201200212 [arXiv:1210.0544 [astro-ph.CO]].[2] N. Aghanim et al. [Planck], Astron. Astrophys. , A6 (2020) doi:10.1051/0004-6361/201833910 [arXiv:1807.06209 [astro-ph.CO]].[3] M. Bauer and T. Plehn, Lect. Notes Phys. , pp. (2019) doi:10.1007/978-3-030-16234-4[arXiv:1705.01987 [hep-ph]].[4] E. Aprile et al. [XENON], Phys. Rev. Lett. , no.11, 111302 (2018)doi:10.1103/PhysRevLett.121.111302 [arXiv:1805.12562 [astro-ph.CO]].[5] X. Ren et al. [PandaX-II], Phys. Rev. Lett. , no.2, 021304 (2018)doi:10.1103/PhysRevLett.121.021304 [arXiv:1802.06912 [hep-ph]].[6] E. Aprile et al. [XENON], Phys. Rev. Lett. , no.14, 141301 (2019)doi:10.1103/PhysRevLett.122.141301 [arXiv:1902.03234 [astro-ph.CO]].[7] S. Knapen, T. Lin and K. M. Zurek, Phys. Rev. D , no.11, 115021 (2017)doi:10.1103/PhysRevD.96.115021 [arXiv:1709.07882 [hep-ph]].[8] T. Lin, PoS , 009 (2019) doi:10.22323/1.333.0009 [arXiv:1904.07915 [hep-ph]].[9] M. Ablikim et al. [BESIII], Nucl. Instrum. Meth. A , 345-399 (2010)doi:10.1016/j.nima.2009.12.050 [arXiv:0911.4960 [physics.ins-det]].[10] B. Aubert et al. [BaBar], Nucl. Instrum. Meth. A , 615-701 (2013)doi:10.1016/j.nima.2013.05.107 [arXiv:1305.3560 [physics.ins-det]].[11] I. Adachi et al. [Belle-II], Nucl. Instrum. Meth. A , 46-59 (2018)doi:10.1016/j.nima.2018.03.068 12] E. Kou et al. [Belle-II], PTEP , no.12, 123C01 (2019) [erratum: PTEP , no.2,029201 (2020)] doi:10.1093/ptep/ptz106 [arXiv:1808.10567 [hep-ex]].[13] J. D. Bjorken, R. Essig, P. Schuster and N. Toro, Phys. Rev. D , 075018 (2009)doi:10.1103/PhysRevD.80.075018 [arXiv:0906.0580 [hep-ph]].[14] A. Boveia and C. Doglioni, Ann. Rev. Nucl. Part. Sci. , 429-459 (2018) doi:10.1146/annurev-nucl-101917-021008 [arXiv:1810.12238 [hep-ex]].[15] M. Pospelov, A. Ritz and M. B. Voloshin, Phys. Lett. B , 53-61 (2008)doi:10.1016/j.physletb.2008.02.052 [arXiv:0711.4866 [hep-ph]].[16] N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer and N. Weiner, Phys. Rev. D , 015014(2009) doi:10.1103/PhysRevD.79.015014 [arXiv:0810.0713 [hep-ph]].[17] M. Pospelov, Phys. Rev. D , 095002 (2009) doi:10.1103/PhysRevD.80.095002[arXiv:0811.1030 [hep-ph]].[18] S. Baek, P. Ko and W. I. Park, JHEP , 013 (2013) doi:10.1007/JHEP07(2013)013[arXiv:1303.4280 [hep-ph]].[19] B. Holdom, Phys. Lett. B , 196-198 (1986) doi:10.1016/0370-2693(86)91377-8[20] E. Izaguirre, G. Krnjaic and B. Shuve, Phys. Rev. D , no.6, 063523 (2016)doi:10.1103/PhysRevD.93.063523 [arXiv:1508.03050 [hep-ph]].[21] M. Duerr, T. Ferber, C. Hearty, F. Kahlhoefer, K. Schmidt-Hoberg and P. Tunney, JHEP ,039 (2020) doi:10.1007/JHEP02(2020)039 [arXiv:1911.03176 [hep-ph]].[22] Y. Zhang, W. T. Zhang, M. Song, X. A. Pan, Z. M. Niu and G. Li, Phys. Rev. D , no.11,115016 (2019) doi:10.1103/PhysRevD.100.115016 [arXiv:1907.07046 [hep-ph]].[23] J. P. Lees et al. [BaBar], Phys. Rev. Lett. , no.13, 131804 (2017)doi:10.1103/PhysRevLett.119.131804 [arXiv:1702.03327 [hep-ex]].[24] G. Mohlabeng, Phys. Rev. D , no.11, 115001 (2019) doi:10.1103/PhysRevD.99.115001[arXiv:1902.05075 [hep-ph]].[25] D. Tucker-Smith and N. Weiner, Phys. Rev. D , 043502 (2001)doi:10.1103/PhysRevD.64.043502 [arXiv:hep-ph/0101138 [hep-ph]].[26] S. Baek, P. Ko and W. I. Park, Phys. Lett. B , 255-259 (2015)doi:10.1016/j.physletb.2015.06.002 [arXiv:1407.6588 [hep-ph]].[27] P. Ko, T. Matsui and Y. L. Tang, JHEP , 082 (2020) doi:10.1007/JHEP10(2020)082[arXiv:1910.04311 [hep-ph]]. 28] S. Baek, J. Kim and P. Ko, Phys. Lett. B , 135848 (2020)doi:10.1016/j.physletb.2020.135848 [arXiv:2006.16876 [hep-ph]].[29] K. Griest and D. Seckel, Phys. Rev. D , 3191-3203 (1991) doi:10.1103/PhysRevD.43.3191[30] T. R. Slatyer, Phys. Rev. D , no.2, 023527 (2016) doi:10.1103/PhysRevD.93.023527[arXiv:1506.03811 [hep-ph]].[31] K. Harigaya, Y. Nakai and M. Suzuki, Phys. Lett. B , 135729 (2020)doi:10.1016/j.physletb.2020.135729 [arXiv:2006.11938 [hep-ph]].[32] N. F. Bell, J. B. Dent, B. Dutta, S. Ghosh, J. Kumar and J. L. Newstead, Phys. Rev. Lett. , no.16, 161803 (2020) doi:10.1103/PhysRevLett.125.161803 [arXiv:2006.12461 [hep-ph]].[33] H. M. Lee, [arXiv:2006.13183 [hep-ph]].[34] M. Baryakhtar, A. Berlin, H. Liu and N. Weiner, [arXiv:2006.13918 [hep-ph]].[35] J. Bramante and N. Song, Phys. Rev. Lett. , no.16, 161805 (2020)doi:10.1103/PhysRevLett.125.161805 [arXiv:2006.14089 [hep-ph]].[36] H. An and D. Yang, [arXiv:2006.15672 [hep-ph]].[37] R. Catena, J. Conrad, C. D¨oring, A. D. Ferella and M. B. Krauss, Phys. Rev. D , no.2,023007 (2018) doi:10.1103/PhysRevD.97.023007 [arXiv:1706.09471 [hep-ph]].[38] H. C. Cheng, Z. Han, I. W. Kim and L. T. Wang, JHEP , 122 (2010)doi:10.1007/JHEP11(2010)122 [arXiv:1008.0405 [hep-ph]].[39] T. Melia, JHEP , 143 (2012) doi:10.1007/JHEP01(2012)143 [arXiv:1110.6185 [hep-ph]].[40] L. Edelhauser, K. T. Matchev and M. Park, JHEP , 006 (2012)doi:10.1007/JHEP11(2012)006 [arXiv:1205.2054 [hep-ph]].[41] N. D. Christensen and D. Salmon, Phys. Rev. D , no.1, 014025 (2014)doi:10.1103/PhysRevD.90.014025 [arXiv:1311.6465 [hep-ph]].[42] S. Y. Choi, Phys. Rev. D , no.11, 115037 (2018) doi:10.1103/PhysRevD.98.115037[arXiv:1811.10377 [hep-ph]].[43] W. Abdallah, A. Hammad, S. Khalil and S. Moretti, Phys. Rev. D , no.9, 095006 (2019)doi:10.1103/PhysRevD.100.095006 [arXiv:1907.08358 [hep-ph]].[44] C. L. Shan, New J. Phys. , 105013 (2009) doi:10.1088/1367-2630/11/10/105013[arXiv:0903.4320 [hep-ph]].[45] B. J. Kavanagh and A. M. Green, Phys. Rev. Lett. , no.3, 031302 (2013)doi:10.1103/PhysRevLett.111.031302 [arXiv:1303.6868 [astro-ph.CO]]. 46] H. C. Cheng and J. Gu, JHEP , 094 (2011) doi:10.1007/JHEP10(2011)094 [arXiv:1109.3471[hep-ph]].[47] L. A. Harland-Lang, C. H. Kom, K. Sakurai and W. J. Stirling, Eur. Phys. J. C , 1969(2012) doi:10.1140/epjc/s10052-012-1969-2 [arXiv:1110.4320 [hep-ph]].[48] L. A. Harland-Lang, C. H. Kom, K. Sakurai and W. J. Stirling, Phys. Rev. Lett. , 181805(2012) doi:10.1103/PhysRevLett.108.181805 [arXiv:1202.0047 [hep-ph]].[49] A. C. Kobach, Phys. Rev. D , 116001 (2013) doi:10.1103/PhysRevD.88.116001[arXiv:1308.5671 [hep-ph]].[50] L. A. Harland-Lang, C. H. Kom, K. Sakurai and M. Tonini, JHEP , 175 (2014)doi:10.1007/JHEP06(2014)175 [arXiv:1312.5720 [hep-ph]].[51] N. D. Christensen, T. Han, Z. Qian, J. Sayre, J. Song and Stefanus, Phys. Rev. D , 114029(2014) doi:10.1103/PhysRevD.90.114029 [arXiv:1404.6258 [hep-ph]].[52] P. Konar and A. K. Swain, Phys. Rev. D , no.1, 015021 (2016)doi:10.1103/PhysRevD.93.015021 [arXiv:1509.00298 [hep-ph]].[53] P. Ko and H. Yokoya, JHEP (2016), 109 doi:10.1007/JHEP08(2016)109 [arXiv:1603.04737[hep-ph]].[54] Q. F. Xiang, X. J. Bi, Q. S. Yan, P. F. Yin and Z. H. Yu, Phys. Rev. D , no.7, 075037(2017) doi:10.1103/PhysRevD.95.075037 [arXiv:1610.03372 [hep-ph]].[55] T. Kamon, P. Ko and J. Li, Eur. Phys. J. C , no.9, 652 (2017) doi:10.1140/epjc/s10052-017-5240-8 [arXiv:1705.02149 [hep-ph]].[56] B. Dutta, T. Kamon, P. Ko and J. Li, Eur. Phys. J. C , no.7, 595 (2018)doi:10.1140/epjc/s10052-018-6071-y [arXiv:1712.05123 [hep-ph]].[57] Z. Flowers, Q. Meier, C. Rogan, D. W. Kang and S. C. Park, JHEP , 132 (2020)doi:10.1007/JHEP03(2020)132 [arXiv:1903.05825 [hep-ph]].[58] D. Kim, K. T. Matchev and P. Shyamsundar, JHEP , 154 (2019)doi:10.1007/JHEP10(2019)154 [arXiv:1906.02821 [hep-ph]].[59] S. Banerjee, B. Bhattacherjee, A. Goudelis, B. Herrmann, D. Sengupta and R. Sengupta,[arXiv:1912.06669 [hep-ph]].[60] K. J. Bae, M. Park and M. Zhang, Phys. Rev. D , no.11, 115036 (2020)doi:10.1103/PhysRevD.101.115036 [arXiv:2001.02142 [hep-ph]].[61] M. Duerr, T. Ferber, C. Garcia-Cely, C. Hearty and K. Schmidt-Hoberg, [arXiv:2012.08595 hep-ph]].[62] C. Boehm, M. J. Dolan and C. McCabe, JCAP , 027 (2012) doi:10.1088/1475-7516/2012/12/027 [arXiv:1207.0497 [astro-ph.CO]].[63] P. F. Depta, M. Hufnagel, K. Schmidt-Hoberg and S. Wild, JCAP , 029 (2019)doi:10.1088/1475-7516/2019/04/029 [arXiv:1901.06944 [hep-ph]].[64] G. F. Giudice, D. Kim, J. C. Park and S. Shin, Phys. Lett. B , 543-552 (2018)doi:10.1016/j.physletb.2018.03.043 [arXiv:1712.07126 [hep-ph]].[65] P. Ilten, Y. Soreq, M. Williams and W. Xue, JHEP , 004 (2018)doi:10.1007/JHEP06(2018)004 [arXiv:1801.04847 [hep-ph]].[66] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer,P. Torrielli and M. Zaro, JHEP , 079 (2014) doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301[hep-ph]].[67] P. A. Zyla et al. [Particle Data Group], PTEP , no.8, 083C01 (2020)doi:10.1093/ptep/ptaa104[68] A. Birkedal, K. Matchev and M. Perelstein, Phys. Rev. D , 077701 (2004)doi:10.1103/PhysRevD.70.077701 [arXiv:hep-ph/0403004 [hep-ph]].[69] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, Comput. Phys. Commun. , 2250-2300 (2014) doi:10.1016/j.cpc.2014.04.012 [arXiv:1310.1921 [hep-ph]].[70] V. V. Gligorov, S. Knapen, M. Papucci and D. J. Robinson, Phys. Rev. D , no.1, 015023(2018) doi:10.1103/PhysRevD.97.015023 [arXiv:1708.09395 [hep-ph]].[71] A. Berlin, S. Gori, P. Schuster and N. Toro, Phys. Rev. D , no.3, 035011 (2018)doi:10.1103/PhysRevD.98.035011 [arXiv:1804.00661 [hep-ph]].[72] J. L. Feng, I. Galon, F. Kling and S. Trojanowski, Phys. Rev. D , no.3, 035001 (2018)doi:10.1103/PhysRevD.97.035001 [arXiv:1708.09389 [hep-ph]].[73] J. P. Chou, D. Curtin and H. J. Lubatti, Phys. Lett. B , 29-36 (2017)doi:10.1016/j.physletb.2017.01.043 [arXiv:1606.06298 [hep-ph]].[74] A. Hook, E. Izaguirre and J. G. Wacker, Adv. High Energy Phys. , 859762 (2011)doi:10.1155/2011/859762 [arXiv:1006.0973 [hep-ph]]., 859762 (2011)doi:10.1155/2011/859762 [arXiv:1006.0973 [hep-ph]].