Exploring properties of long-lived particles in inelastic dark matter models at Belle II

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, the co-annihilation becomes the dominant channel for thermal relic density and the DM excited state can be long-lived at the collider scale. We study scalar and fermion inelastic dark matter models for O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(1) GeV DM at Belle II with U(1)D dark gauge symmetry broken into its Z2 subgroup. We focus on dilepton displaced vertex signatures from decays of the DM excited state. With the help of precise displaced vertex detection ability at Belle II, we can explore the DM spin, mass and mass splitting between DM excited and ground states. Especially, we show scalar and fermion DM candidates can be discriminated and the mass and mass splitting of DM sector can be determined within the percentage of deviation for some benchmark points. Furthermore, the allowed parameter space to explain the excess of muon (g− 2)μ is also studied and it can be covered in our displaced vertex analysis during the early stage of Belle II experiment.


Introduction
The nature of dark matter (DM) in our Universe is still a great mysterious issue. As we know, DM plays a major role in the structure formation [1], and its abundance is about 5.5 times larger than the ordinary matter in the present Universe [2]. However, even all robust evidences to support the existence of DM until now are only connected with gravitational interactions. 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 Interacting Massive Particle (WIMP) DM candidate with the freezeout mechanism has been overwhelming in both theoretical and experimental communities during the past decades as shown in ref. [3] and references therein. This kind of WIMP DM particles can be searched with direct, indirect detections and also at colliders experiments.

JHEP04(2021)269
The null signal result from direct detection provides severe bounds on the cross section of DM and nuclear scattering above a few GeV [4][5][6]. For instance, the spin independent (spin dependent) DM and nuclear scattering cross section is restricted to σ SI 10 −46 cm 2 [4] (σ SD 2 × 10 −41 cm 2 [6]) for M DM = 100 GeV. Nevertheless, lighter DM candidates are still less constrained because the restrictions of fine energy threshold are required from these direct detection experiments. Therefore, those dark sector models with MeV to GeV DM candidates become more and more 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 powerful than the high-energy machines such as the Tevatron and LHC [14]. On the other hand, DM is not the lonely particle in the dark sector for most of DM portal models [15][16][17][18]. In order to connect the SM sector with the dark sector, a mediator is required. Among these mediator candidates, the dark photon via the kinetic mixing portal is an attractive type [19] for the dark sector with Abelian gauge symmetry. In order to be consistent with the relic density bound, 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 constraints can be set up as shown in refs. [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 them to explore some models with light DM candidates via the monophoton signature [20][21][22][23]. Furthermore, as pointed out in refs. [20,21,24], we can not only study the mono-photon signature, 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 the invisible DM signature cannot reach at high-intensity machines. Hence, we will focus on exploring the properties of long-lived particles in dark sectors with U(1) D gauge symmetry at Belle II in this work.
The inelastic (or excited) DM models with extra U(1) D gauge symmetry [25][26][27][28] is one of the most popular dark sector models with light DM candidates. 1 There are at least two states 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 dominant one of DM relic density in early Universe [29]. It is one of the unique features of this sort of models. The co-annihilation production for light DM 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 constraint from DM and nuclear inelastic scattering is much weaker than the elastic one in the direct detection experiments. 2 It makes more allowed parameter space can be explored in the inelastic DM models. Besides, there are also other rich phenomenons 1 In ref. [25], the U(1)D symmetry is explicitly and softly broken by a dim-2 operator for scalar DM, and there is no dark Higgs boson there. On the other hand, in the models considered in this paper in that there is dark Higgs boson which plays important roles in DM phenomenology. 2 The bounds from direct detection experiments are weak for ∼ O(1) GeV DM. Hence we can safely ignore possible contributions to the DM-nucleon elastic scattering from the SM Higgs boson and dark Higgs boson in the direct detection experiments.

JHEP04(2021)269
in this model. For example, it is possible to explain the muon g − 2 anomaly [24] and XENON1T excess [28,[31][32][33][34][35][36][37]. Last but not least, the excited DM state can naturally become long-lived and leave displaced vertex inside detectors after it has been produced such that these novel signatures can be searched for at colliders. The fundamental properties of DM particle include its mass, charge, and spin. Since we assume DM is electrically neutral or millicharged, only its mass and spin are left to be determined. Actually, the mass range for WIMP DM can cover from few MeV to a few hundred TeV and its spin can be 0, 1/2, 1, 3/2 or even higher. Most of the time the precise determination of DM spin and mass are challenging tasks. However, there are already some previous studies about DM spin and mass determinations. For the DM spin determination, it has been studied for direct detection experiments in ref. [38] and collider experiments in refs. [39][40][41][42][43][44]. For the DM mass determination, it has been studied for direct detection in refs. [45,46] and collider experiments in refs. [47][48][49][50][51][52][53][54][55][56][57][58][59][60][61]. Especially, using the displaced vertex information for DM mass determination starts from ref. [58]. Thanks to the precise track resolution in the inner detector of ATLAS and CMS experiments, the DM mass reconstruction can be studied in some BSM models of the specific cascade processes involving long-lived particles. Furthermore, the resolutions of displaced vertex, charged lepton and photon momentum can be even improved at detectors of Belle II experiments [11]. Therefore, our goal in this paper is to explore the spin and mass properties of inelastic DM models with dilepton displaced vertex signatures at Belle II. As we will see in section 4, the scalar and fermion inelastic DM models can be well discriminated and the DM mass and mass splitting between DM excited and ground states can be determined within the percentage of deviation for those benchmark points.
The organization of this paper is as follows. We first review scalar and fermion inelastic DM models with U(1) D gauge symmetry in section 2. Both analytical representations and numerical results for cross sections of relevant signal processes are displayed in section 3. We also point out how to distinguish scalar and fermion inelastic DM models via the size and kinematic distribution of their cross sections in this section. Detailed simulations and methods to determine DM mass, and mass splitting between DM excited and ground states are shown in section 4. Finally, we conclude our studies in section 5. Some supplemental formulae for section 3 can be found in the appendix A and B.

Inelastic dark matter models
The scalar and fermion inelastic (or excited) DM models with U(1) D gauge symmetry are reviewed in this section. After the spontaneous symmetry breaking (SSB) of this U(1) D gauge symmetry, we expect the accidentally residual Z 2 symmetry, φ 1 → −φ 1 (scalar) or χ 1 → −χ 1 (fermion), can be left such that φ 1 or χ 1 are stable and become DM candidates in our Universe.

JHEP04(2021)269
U(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 particles do not carry U(1) D charges.
The scalar part of the renormalizable and gauge invariant Lagrangian density is where W a µ , B µ , and X µ are the gauge potentials of the SU(2) L , U(1) Y and U(1) D with gauge couplings g, g and g D , respectively. The σ a is the Pauli matrix and a runs from 1 to 3. The scalar potential in eq. (2.1) is given by where all parameters are assumed to be real for simplicity. We then expand H, Φ fields around the vacuum with the unitary gauge, and highlight some important relations in the below: • The four-points and off-diagonal interactions of the DM sector and new gauge boson can be obtained from the |D µ φ| 2 term: • The DM sector mass splitting and dark Higgs, DM sector trilinear interactions are derived from µ Φφ Φ * φ 2 + H.c. terms: • Finally, the φ 1 , φ 2 masses and their mass splitting can be represented as and

JHEP04(2021)269
Notice the interaction states (h, h D ) would be rotated to the mass states (h 1 , h 2 ) via the mixing angle θ. The sin θ is chosen to be small enough in our analysis such that it is consistent with the SM-like Higgs boson data at the LHC. 3 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 µν . The Lagrangian density of this part can be represented as where is the kinetic mixing parameter between these two U(1)s. If we apply the linear order approximation in , the extra interaction terms for SM fermions and Z boson can be written as where c W is the weak mixing angle and x l = −1, x ν = 0, x q = 2 3 or − 1 3 depending on the electrical charge of quark. The Z boson mass can be approximated as Notice the correction from the kinetic mixing term is second order in which can be safely neglected here.

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 Q D (Φ) = +2 and Q D (χ) = +1. Again, the SM-like Higgs doublet and other SM particles do not carry U(1) D charges. The scalar part of the renormalizable and gauge invariant Lagrangian density is where D µ H and D µ Φ are the same in eq. (2.2) and (2.3). The scalar potential in eq. (2.14) is given by where all parameters are assumed to be real for simplicity. The Lagrangian density of dark matter sector part is -5 -

JHEP04(2021)269
where f is assumed to be a real parameter. In order to decompose the Dirac fermion χ into a pair of two independent Majorana fermions, χ 1 and χ 2 , we set After expanding H, Φ fields around the vacuum with the unitary gauge as shown in eq. (2.6), the eq. (2.16) can be written as where χ 1 , χ 2 masses and their mass splitting can be represented as and Finally, the Z boson mass and its interactions with SM fermions are the same with eq. (2.13) and (2.12).

The relevant cross sections
In this section, we show both analytical representations and numerical results for cross sections of e + e − → φ 1 φ 2 (χ 1 χ 2 ) and e + e − → φ 1 φ 2 (χ 1 χ 2 )γ processes. Especially, the method to distinguish spin-0 and spin-1/2 DM candidates in the inelastic DM models are also discussed here.

The analytical representations
The differential cross section for e + e − → φ 1 φ 2 via the s-channel Z can be represented as where s, t are Mandelstam variables and Γ Z is the total width of Z boson. In order to study the differential angular distributions of cross sections, we transfer from dσ/dt to dσ/d cos θ, where θ is the polar angle of φ 2 and it is defined as the direction of φ 2 relative to the positron beam direction.
The dσ/d cos θ in the centre-of-mass (CM) frame can be simply written as and where E CM is the centre-of-mass energy. We can find eq. (3.2) is always proportional to 1 − cos 2 θ. Consequently, it indicates if the s-channel Z is on-shell produced, it is longitudinally polarized (helicity = 0). However, the formula for dσ/d cos θ in the Belle II laboratory (LAB) frame is more tedious, so we show it in eq. (A.3) 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 to the initial electron beam direction for the Belle II machine compared with the ones in the CM frame.
The differential cross section for e + e − → χ 1 χ 2 via s-channel Z can be represented as and dσ/d cos θ in the CM frame can be written as where ξ can be found in eq. (3.3) with M φ 1,2 → M χ 1,2 . For the M χ 1,2 → 0 limit, eq. (3.5) is proportional to 1 + cos 2 θ, which is well known. It shows if the s-channel Z is on-shell produced, it is transversely polarized (helicity = ±). However, the mass effects of M χ 1,2 spoil parts of this property. Similarly, we show the formula for dσ/d cos θ in the LAB frame in eq. (A.4) of the appendix A. Again, there is a skewing behavior for the differential angular distributions 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 energy or soft objects at lepton colliders. Especially, for the process with on-shell Z production, the ISR photon is used not only for background rejection, but also for Z invariant mass reconstruction. To include the ISR photon in e + e − → φ 1 φ 2 (χ 1 χ 2 ) process, the differential cross section can be written as [69] dσ and the splitting kernal P(z, cos θ ) is where α is the fine structure constant, z = Eγ E CM /2 is the energy fraction of ISR photon from the initial electron/positron and θ is the polar angle of ISR photon and it is defined as the direction of γ to the positron beam direction. The differential form of σ(e + e − → φ 1 φ 2 (χ 1 χ 2 )) can be found in eq. (3.1) and (3.4). In the CM frame, if we assign the four momentum of the initial electron/positron and ISR photon as p e ± = (E CM /2, 0, 0, ±E CM /2) and p ISR = (E γ , p x,γ , p y,γ , p z,γ ), the four momentum of electron/positron after the radiation is . Let's consider two limit scenarios to intuitively catch up the behavior of d σ(e + e − → φ 1 φ 2 (χ 1 χ 2 ))/d cos θ distributons. First, if the ISR photon is soft (z ≈ 0), p e ± ≈ p e ± and d σ(e + e − → φ 1 φ 2 (χ 1 χ 2 ))/d cos θ distributions are close to eq. (3.2) and (3.5). On the contrary, if the ISR photon takes almost all energy from the initial electron/positron (z → 1), the φ 2 (χ 2 ) is obvious in the forward or backward direction. Since d σ(e + e − → φ 1 φ 2 (χ 1 χ 2 ))/d cos θ distributions are highly dependent on the z parameter of the ISR photon, the e + e − → φ 1 φ 2 (χ 1 χ 2 )γ processes are not ideal to distinguish scalar and fermion inelastic DM models.
Finally, we focus on the parameter space

The numerical results
We first generate both scalar and fermion inelastic DM UFO model files from Feyn-Rules [70], and then calculate cross sections of e + e − → φ 2 φ 1 (χ 2 χ 1 ) and e + e − → φ 2 φ 1 (χ 2 χ 1 )γ processes via MadGraph5 aMC@NLO. 6 The beam energies are set to be E(e + ) = 4.0 GeV and E(e − ) = 7.0 GeV which are consistent with the Belle II experi- 6 The Z boson total decay width is automatically calculated in MadGraph5 aMC@NLO.
ment. In order to display the relations of cross sections with M φ 1 ,χ 1 and m Z , we study processes are shown in the left panel of figure 1. Since all of these cross sections are proportional to 2 α D , it is straightforward to rescale cross sections with different values of and α D . On the other hand, the influence from the changes of ∆ φ,χ to cross sections are mild for ∆ φ,χ < 0.5M φ 1 ,χ 1 . The scalar and fermion pair production cross sections can be scaled by β 3/2 and β 1/2 respectively, where β is the velocity of the final state particle in the CM frame. Because of this extra β factor for the scalar case, cross sections for e + e − → φ 2 φ 1 are suppressed compared with e + e − → χ 2 χ 1 .
On the other hand, the relation between M φ 1 ,χ 1 (GeV) and σ (pb) for e + e − → φ 2 φ 1 (χ 2 χ 1 )γ processes are shown in the right panel of figure 1 with the same parameter settings. Here, the basic cuts E(γ) > 0.1 GeV and η < 2.203 are applied for the ISR photon. We find the main contributions of e + e − → φ 2 φ 1 (χ 2 χ 1 )γ come from e + e − → Z γ → φ 2 φ 1 (χ 2 χ 1 )γ processes. It explains why scalar and fermion pair cross sections in this process are very close to each other in the right panel of figure 1.
We then turn to the study of (1/σ)(dσ/d cos θ) distributions for e + e − → φ 2 φ 1 and e + e − → χ 2 χ 1 processes with fixed α D = 0.1, = 0.01, ∆ φ,χ = 0.1M φ 1 ,χ 1 , and m Z = 3M φ 1 ,χ 1 . The results in the CM and Belle II LAB frame are shown in figure 2. It is clear that the (1/σ)(dσ/d cos θ) distributions for e + e − → φ 2 φ 1 in the CM frame is proportional to 1 − cos 2 θ and e + e − → χ 2 χ 1 behaves close to 1 + cos 2 θ with some distortions from M χ 1,2 effects. However, since the measurement of time-of-flight for the long-lived particle is poor at the Belle II machine, we will lose this information and cannot make the Lorentz transformation for 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 in the Belle II LAB frame. Hopefully, the real situation for the distributions in the Belle II LAB frame is just shifted to the electron beam direction for the Belle II machine compared with the ones in the CM frame. It is still clear to see the different distributions between scalar and fermion inelastic DM models.  Finally, we fix the same parameter settings as figure 2 for the (1/σ)(dσ/d cos θ) distributions of e + e − → φ 2 φ 1 (χ 2 χ 1 )γ processes in the CM and LAB frames in figure 3. Again, the distributions in the CM frame is only for the comparison and it is clear to see the skewing behavior for the differential angular distributions of cross sections in the LAB frame compared with the ones in the CM frame. As discussed in the previous subsection, because these distributions are highly dependent on the z parameter of the ISR photon, the e + e − → φ 1 φ 2 (χ 1 χ 2 )γ processes are not ideal for the DM spin discrimination. Only when the Z is on-shell produced and heavy, the ISR photon becomes soft (z is small) such that differences of (1/σ)(dσ/d cos θ) distributions from scalar and fermion inelastic DM models show up.
In summary, the size for cross sections of e + e − → φ 2 φ 1 (χ 2 χ 1 ) and their (1/σ)(dσ/d cos θ) distributions in the Belle II LAB frame can help us to statistically distinguish fermion inelastic DM model from the scalar one even for the same parameter settings. and then show some interesting kinematic distributions based on four signal benchmark points. We further set up event selections for dilepton displaced signature and relevant results are discussed. Finally, we solve kinematic equations of e + e − → χ 1 χ 2 → χ 1 χ 1 l + l − and e + e − → Z γ → χ 1 χ 2 γ → χ 1 χ 1 l + l − γ processes to determine both M χ 1 and M χ 2 .

The Belle II experiment and signal signatures
The SuperKEKB accelerator of Belle II experiment is a circular asymmetric e + e − collider with the nominal collision energy of √ s = 10.58 GeV. The beam parameters are E(e + ) = 4 GeV and E(e − ) = 7 GeV. The planned full integrated luminosity for the final dataset is 50 ab −1 . In this work, the following Belle II sub-detectors are relevant: the tracking system including vertex detectors (VXD) and central drift chamber (CDC), the electromagnetic calorimeter (ECAL), and muon system.
We consider the following two kinds of processes in inelastic DM models: 7 Because the major contribuitons of e + e − → φ1φ2(χ1χ2)γ come from the on-shell Z production, hence, we only focus on the process e + e − → Z γ → φ1φ2(χ1χ2)γ in the following analysis unless noted otherwise.

JHEP04(2021)269
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 products are too soft to be detected in eq. (4.2).
• Mono-γ with prompt lepton pair: the excited DM state promptly decays and the visible products can be successfully detected and defined in eq. (4.2).
• Mono-γ with displaced lepton pair: the excited DM state is long-lived and decays inside the detector leaving the displaced vertex in eq. (4.2).
• Only prompt lepton pair: the same as the second one but without ISR photon in eq. (4.1).
• Only displaced lepton pair: the same as the third one but without ISR photon in eq. (4.1).
The Mono-γ signature is well-studied as shown in refs. [21][22][23]. Since we focus on searching for long-lived particles in this paper, only signatures with displaced vertices are studied.
On the other hand, the signature with displaced charged pion pairs is also possible from the long-lived excited DM state decay. However, for simplicity, here we only study signatures with displaced electron and muon pairs. In order to make our results more realistic, we follow ref. [11] to involve the detector resolution effects with Gaussian smearing at Belle II experiment. The tracking resolution of electron/muon momentum in the CDC is given by where p l t is the transverse momentum of electron/muon track and β is its velocity in the natural unit. We conservatively apply σ p l /p l = 0.005 in our event analysis. On the other hand, the muon efficiency in the muon system is approximated to 0.98. The photon momentum resolution in the ECAL is approximated to σ Eγ /E γ = 2% where E γ is the energy of photon. Finally, we use the resolution of σ r DV = 26µm for the displaced vertex vector of lepton pair in our analysis.

Benchmark points and kinematic distributions
In this study, we use the following four benchmark points (BPs) to display our analysis:  with fixed α D = 0.1. 8 The first BP is inspired from the allowed parameter space to explain the muon anomalous magnetic moment excess [24]. The other three BPs are the ones which cannot be covered from Mono-γ searches from BaBar and Belle II [21]. In addition, they can also be searched for in the future long-lived particle experiments, like CODEX-b [71], SeaQuest [72], FASER [73], MATHUSLA [74], SHiP [75], ANUBIS [76] and AL3X [77].
We show some interesting kinematic distributions based on the above four BPs for e + e − → χ 1 χ 2 → χ 1 χ 1 l + l − and e + e − → Z γ → χ 1 χ 2 γ → χ 1 χ 1 l + l − γ in figure 4 and figure 5, respectively. First, the lepton pair invariant mass is smaller or equal to the mass splitting ∆ χ (dashed lines), 2m l ≤ M ll ≤ ∆ χ . Once we have enough signal events, the threshold value of M ll can help us to roughly determine the mass splitting ∆ χ in inelastic DM models. Second, the distance of displaced vertex of χ 2 is not only dependent on M χ 2 , ∆ χ , and α D , but also the boost of χ 2 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 χ 2 decay length in BP1 (BP2) is the shortest (longest) one on average. The distributions for projection of the decay vertex distance on the transverse plane, R xy , are similar. On the other hand, E ll and ∆R ll are energy and angular distance for the lepton pair. We can find the boost of χ 2 is smaller for the process involving ISR photon by comparing r z , E ll and ∆R ll distributions in figure 4 8 Notice Z → φ1φ2 or χ1χ2 is the dominant decay channel for these four BPs since Z can only couple to SM fermion pairs via the kinetic mixing. We have checked the branching ratio for Z to SM particles is less than 10 −3 for these four BPs. Therefore, we can safely ignore relevant constraints from visible dark photon searches.  and figure 5. The larger boost on χ 2 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 detector resolutions such that it is more challenge to reconstruct DM mass and mass splitting between DM excited and ground states in the real experiment. Third, if the Z is light enough, the ISR photon and Z are almost back-to-back produced such that the lepton from Z decay can have large angular distance with the ISR photon. Here we show the maximum angular distance between ISR photon and one of the lepton in lepton pair, ∆R max γl . Fourth, we boost the E(γ) distribution of ISR photon from the Belle II LAB frame to the CM frame and it is clear to see its monoenergetic behavior. This feature can largely reduce SM backgrounds which are continuous in the photon energy spectrum. Finally, the Z invariant mass can be reconstructed from  eq. (D.6) with measurable observables. Because of the detector resolution, it is clear to find smearing effects from the true Z mass (dashed lines), especially for the case of light Z . We will see this phenomenon critically affects abilities to determine DM mass and the mass splitting between DM excited and ground states for the process in eq. (4.2). Notice those kinematic distributions for the scalar inelastic DM model are similar to figure 4 and figure 5, so we do not show them again.

Event selections and results
We closely follow ref. [21] to set up event selections for the displaced signature at Belle II. We only conservatively consider the following two background-free regions after event selections in our analysis: • low R xy region (100% detection efficiency): 0.2 < R xy ≤ 0.9 cm (electron), 0.2 < R xy ≤ 17.0 cm (muon).
According to the arguments in ref. [21], backgrounds from photon conversion can be further reduced to a negligible level by requiring the lepton pair invariant mass M ll ≥ 0.03 GeV, and their opening angle larger than 0.1 rad. Notice backgrounds from the photon conversion to muon pair in the region 0.9 < R xy ≤ 17.0 cm are negligible such that the analysis for muon in low R xy region can be extended. We summarize all event selections in our analysis in table 1. For the muon pair in the final state, we veto the invariant mass region, 0.48 ≤ m µ + µ − ≤ 0.52 GeV, to reject backgrounds from K 0 S decay. More discussions on the trigger issue can be found in ref. [21]. Here, event selections in table 1 are done offline and we simply assume these events are already triggered and stored.
The results for four BPs are shown in table 2 and 3 for processes in eq. (4.1) and (4.2), 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 −1 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. The discrimination of scalar and fermion inelastic DM models with N event in table 2 and 3 is shown in appendix C.
Finally, we use the most optimistic value of 50 ab −1 for event selections in table 1 at Belle II for scalar and fermion inelastic DM models and predict the future bounds from e + e − → φ 1 φ 2 (χ 1 χ 2 ) and e + e − → φ 1 φ 2 (χ 1 χ 2 )γ processes in figure 6. Notice the off-shell Z production in the second process is also considered. Here we fix the parameters, α D = 0.1, m Z = 3M φ 1 ,χ 1 and ∆ φ,χ = 0.1M φ 1 ,χ 1 , and apply 90% C.L. contours which correspond to an upper limit of 2.3 events with the assumption of background-free. Notice constraints from model-independent LEP bound [78] and BaBar mono-γ bound [23] are added for the comparison. We closely follow ref. [21] for the recasting of BaBar mono-γ constraints in inelastic DM models. The correct relic abundance lines are also shown in figure 6. Here the dominant contribution comes from φ 1 φ 2 (χ 1 χ 2 ) coannihilation process [20]. In our scenario, φ 1 φ 1 (χ 1 χ 1 ) annihilation processes with dark Higgs boson or SM-like Higgs boson are either ignorable or kinematically forbidden. For M φ 1 ,χ 1 < 0.4 GeV, even dilepton displaced vertex is located in our target regions, the lepton pair in the final state is too soft and collinear to pass the event selections in table 1. On the contrary, for M φ 1 ,χ 1 4.9 GeV, the lepton pair in the final state becomes more energetic and well-separated, but dilepton displaced vertex is too short to reach our target regions and production cross sections are also highly suppressed. We find our results in the second process for fermion inelastic DM model are consistent with the ones in ref. [21] and we further include detector resolution effects in our analysis. Besides, we believe our results for the first process in both scalar and fermion -16 -  Table 2. The analysis results of four BPs for the process in eq. (4.1). The first column is the type of inelastic DM models. The second column is labeled by BP which e and µ denote the type of lepton pair in the final state. The third column is production cross sections for process in eq. (4.1). 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 events with an integrated luminosity of 1ab −1 .  e + e -→χ1χ2γ Figure 7. The same as figure 6 but for m Z = 3M φ1,χ1 and ∆ φ,χ = 0.1M φ1,χ1 . The green shaded region bounded by the green dashed lines is the 2σ allowed region for the (g − 2) µ excess and the lighter gray region excluded by the (g − 2) µ at 5σ C.L. inelastic DM models are first shown in the literature. Most importantly, the future bounds from the first process can be even stronger than the usual ones in the second process.

JHEP04(2021)269
Moreover, inspired by BP1 which can explain the muon (g − 2) µ [24], we perform the same analysis in figure 7 for α D = 0.1, m Z = 3M φ 1 ,χ 1 and ∆ φ,χ = 0.4M φ 1 ,χ 1 fixed. We show the 2σ allowed and the 5σ excluded regions for the muon (g −2) µ as well as the modelindependent LEP bound [78], BaBar mono-γ bound [23] and also correct relic abundance lines. Once the mass of DM is increasing, the mass splitting between DM ground and excited states is also enhanced, hence, the final state lepton pair becomes more energetic. Therefore, one can find the BaBar mono-γ bound is slightly weakened when the displaced channel is open. Since the recasting in ref. [24] does not include the requirement on the angle θ DV LAB in table 1, the results are weaker than the ones 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-γ bound and it is shown in dotted gray line in figure 7 which can already cover the (g −2) µ allowed window in the parameter space. Even the recasting BaBar mono-γ bound can close this region, we find our result can give much stronger bound, ∼ O(10 −5 − 10 −4 ), especially in the parameter space where the mono-γ searches are suppressed. Therefore we can explicitly cover this area using displaced lepton search at the Belle II in a very first stage of the run.
In summary, we found the future bounds from the process without ISR photon can be stronger than the one with ISR photon in figure 6. Instead, in the lower M φ 1 region, the bound for process with ISR photon in figure 6 and figure 7 get stronger because its cross sections are larger than the ones without ISR photon as shown in figure 1.

Determination of the DM mass and mass splitting of DM sector
In general, we cannot uniquely determine the DM mass at colliders because of lacking enough constraints for invisible particles in the final state. We can take a mono-γ search -18 -

JHEP04(2021)269
for DM pair production at B-factories as an example. If the DM is denoted by χ, the process for mono-γ search at B-factories is e + e − → χχγ (or e + e − → χχγ). There are eight unknown values from four-momentum of two DMs in the final state. Unfortunately, only five constraints in this process: four from the four-momentum conservation and one from the same mass for the DM pair. We still need three extra conditions to uniquely determine the DM mass for each event.
Now we turn our attention to the case of inelastic DM models. For the process in eq. (4.1) or (4.2), if the φ 2 (χ 2 ) 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 φ 1 s(χ 1 s) in the final state. However, because of the charge neutrality of the φ 2 (χ 2 ), a three-momentum vector of φ 2 (χ 2 ) is proportional to the direction of displaced vertex (DV) [58,61] where − → p φ 2 (χ 2 ) is the three-momentum vector of φ 2 (χ 2 ) and r DV is the unit vector of displaced vertex from φ 2 (χ 2 ). Therefore, we have two more constraints, and there are seven constraints for this kind of processes in total. We still need one more condition to uniquely determine the DM mass. We first show some key kinematic equations for processes in eq. (4.1) and (4.2) and solve them event-by-event from our Monte Carlo samples and then involve detector resolution effects. 9 More detalied derivations of these kinematic equations can be found in appendix D.
For e + e − → χ 1 χ 2 → χ 1 χ 1 l + l − , we first solve the energy of χ 2 for the subprocess e + e − → χ 1 χ 2 as where E − , E + , and θ are e − , e + beam energies and the polar angle of DV. We then turn to the subprocess χ 2 → χ 1 l + l − . With the help of energy and momentum conservation, we can receive the following equation for inputs of E V , − → p V and r DV , where E χ 2 is shown in eq. (4.5). In figure 8, we display three arbitrary Monte Carlo events on the (M χ 1 , ∆ χ ) plane for four BPs according to eq. (4.6). On the other hand, the black line shows the kinematic endpoint measurement M max l + l − . We can find that the lines from all events and M max l + l − cross to the same point which is the true (M χ 1 , ∆ χ ) in our four BPs. This is a simple application of the Kinematic Focus Point Method proposed in ref. [59]. Similarly, for e + e − → Z γ → χ 1 χ 2 γ → χ 1 χ 1 l + l − γ, we can also solve the energy of χ 2 for the subprocess e + e − → Z γ → χ 1 χ 2 γ as Finally, the kinematic equation for χ 2 → χ 1 f f is the same as eq. (4.6) with E χ 2 in eq. (4.7). Again, we show three arbitrary Monte Carlo events on the (M χ 1 , ∆ χ ) plane for four BPs according to eq. (4.6) in figure 9. All events and M max l + l − cross to the true (M χ 1 , ∆ χ ) in our four BPs.
In order to make the Kinematic Focus Point Method fit to reality, we involve the detector resolution effects from section 4.1 and event selections in table. 1. Here we conservatively assume there are 100 signal events for each BP can be recorded. Since we can solve the above kinematic equations for any two signal events, we will get C 100 2 = 4950 solutions from 100 signal events. After removing the unphysical solutions, we show the distributions of these solutions on the (M χ 1 , M χ 2 ) plane for four BPs of e + e − → χ 1 χ 2 → χ 1 χ 1 l + l − and e + e − → Z γ → χ 1 χ 2 γ → χ 1 χ 1 l + l − γ processes in figures 10 and 11, respectively. Here the bin size for both x and y axes are set to be 10 MeV. Moreover, the most probable   In figure 10, except for BP1 on top left panel, M χ 1 and M χ 2 can be determined within the percentage of deviation for other three BPs. As shown in figure 4, both the decay length of χ 2 is relative short and E ll is relative large for BP1 compared with other three BPs such that events in BP1 cannot be avoided to severely suffer from detector resolution effects. This explains why the physical solutions are reduced and the reconstructed M χ 1 and M χ 2 are shifted to larger values with higher smearing for BP1. In figure 11, M χ 1 and M χ 2 can still be determined within the percentage of deviation for BP2 and BP4. However, apart from detector resolution effects from charged leptons and displaced vertex, the ability to precisely pin down M χ 1 and M χ 2 for the process in eq. (4.2) also relies on how well the on-shell Z can be reconstructed. As shown in figure 5, the peak of m Z is rather spread out (or even shifted) for BP1 and BP3 than for BP2 and BP4. According to eq. (D.6), the reconstructed Z 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 the on-shell Z becomes poor. In the end, the physical  solutions of BP1 and BP3 are largely reduced and shifting behaviors of M χ 1 and M χ 2 become severe in figure 11 compared with the ones in figure 10. This is another reason to encourage our experimental colleagues to search for not only the usual process in eq. (4.2) at B-factories, but also the process in eq. (4.1) with the displaced vertex trigger which is more sensitive to determine the mass and mass splitting of DM sector for some parameter space in the inelastic DM models.

Conclusion
The dark sectors with light DM candidate have received substantial attention as the null signal result from DM direct detetion has been reported until now. On the other hand, the null results from beyond the SM searches at the LHC also hint new physics signatures may hide in the elusive corner. New kind of signatures such as long-lived particles at colliders become more and more popular and may guide a royal road for new physics evidences. Among these dark sector models, the inelastic DM model is an appealing example which can both allow light 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 in this work. . The same as figure 10, but the solutions for eq. (4.6) with E χ2 in eq. (4.7) of e + e − → Z γ → χ 1 χ 2 γ → χ 1 χ 1 l + l − γ process.
In the inelastic DM models, if the mass splitting between DM excited and ground states is small enough, the co-annihilation becomes the dominant channel for thermal relic density and the DM excited state can be long-lived at the collider scale. We first review scalar and fermion inelastic DM models with U(1) D gauge symmetry in section 2 and point out the dark Higgs sector caused the mass splitting between DM excited and ground states and provided a mechanism to the Z mass. Besides, the off-diagonal interaction between Z and DM sector is derived from the covariant derivative term of DM sector.
The analytical representations and numerical results for cross sections of e + e − → φ 1 φ 2 (χ 1 χ 2 ) and e + e − → φ 1 φ 2 (χ 1 χ 2 )γ processes are studied in section 3. In the first process, there is an extra β factor suppression for boson pair cross sections compared with fermion pair ones. In the second process, the dominant channel comes from e + e − → Z γ → φ 1 φ 2 (χ 1 χ 2 )γ. Therefore, cross sections in this process for scalar and fermion inelatic DM models are very close to each other. In addition, the polar angle distributions of φ 2 (χ 2 ) at the Belle II LAB frame as shown in figure 2 are helpful to distinguish these two models.
The novel dilepton displaced vertex signatures at Belle II are the main targets of this study. We include detector resolution effects in section 4.1 and event selections in  Table 5. The same as table. 4, but in e + e − → Z γ → χ 1 χ 2 γ → χ 1 χ 1 l + l − γ process. and 7, which indicate that the future bounds from dilepton displaced vertex searches in inelastic DM models at Belle II are stronger than previous constraints for two benchmark mass windows 0.4 M φ 1 ,χ 1 4.9 GeV and 0.1 M φ 1 ,χ 1 1.1 GeV, respectively. For the latter case, our results can cover the allowed parameter space which can explain the muon (g − 2) µ . Therefore the early stage of the Belle II experiment can explicitly close off this area in the parameter space. We further apply the Kinematic Focus Point method for e + e − → χ 1 χ 2 → χ 1 χ 1 l + l − and e + e − → Z γ → χ 1 χ 2 γ → χ 1 χ 1 l + l − γ processes to determine both M χ 1 and M χ 2 masses. Especially, we find the mass and mass splitting of DM sector can be determined within the percentage of deviation as shown in table. 4 and 5. Before closing, we would like to mention that our analysis in this work is quite general and can be applied to other models such as [79][80][81][82][83].