Exposing Dark Sector with Future Z-Factories

We investigate the prospects of searching dark sector models via exotic Z-boson decay at future $e^+ e^-$ colliders with Giga Z and Tera Z options. Four general categories of dark sector models: Higgs portal dark matter, vector portal dark matter, inelastic dark matter and axion-like particles, are considered. Focusing on channels motivated by the dark sector models, we carry out a model independent study of the sensitivities of Z-factories in probing exotic decays. The limits on branching ratios of the exotic Z decay are typically $\mathcal{O} (10^{-6} - 10^{-8.5}) $ for the Giga Z and $\mathcal{O} (10^{-7.5} - 10^{-11})$ for the Tera Z, and they are compared with the projection for the high luminosity LHC. We demonstrate that future Z-factories can provide its unique and leading sensitivity, and highlight the complementarity with other experiments, including the indirect and direct dark matter search limits, and the existing collider limits. Future Z factories will play a leading role to uncover the hidden sector of the universe in the future.

Searching for dark sector particles, including dark matter (DM) itself and other associated states, is a central goal of many experimental programs around the world. In the mass range between MeV and TeV, collider search remains a crucial method to look for these hidden particles. Since the dark sector particles typically only have weak couplings with the Standard Model, colliders with higher luminosity are natural places to lead this quest. Recently, there have been a couple of proposals for future Z-factories based on circular e + e − colliders, including FCC-ee and CEPC [1][2][3][4], which are considering both Giga-Z and Tera-Z options. Giga-Z (Tera-Z) means running the electron collider at Z pole energy and accumulate 10 9 (10 12 ) Z's respectively. Given the measured cross-section of hadronic Z is 30.5 nb [5], the integrated luminosity for Giga Z (10 9 Z) and Tera Z (10 12 Z in the plan of FCC-ee) are 22.9 fb −1 and 22.9 ab −1 , respectively. In this paper, we give projections on the sensitivities of Z-factory searches to a set of Z rare decay channels inspired by the dark sector models.
In addition to using exotic decay measurements to probe these models, we also compare the reach with direct and indirect dark matter detection experiments, current limits from collider searches, and estimated sensitivities of high luminosity run of the LHC (HL-LHC). Our results demonstrate that the Z-factory measurement will provide the leading sensitivities in many cases. We also include thermal relic abundance, with the understanding that it should serve as an interesting benchmark point, rather than a strict limit.
In section II, we briefly outline the DM indirect, direct searches and the DM relic abundance, in order to compare with the Z decay searches in section III. Section III focuses on well-defined and representative dark sector models to illustrate the power of exotic Z decay search at Z-factory. Certainly, we can not cover all the dark sector models for exotic Z decay. Therefore, we list the possible topologies for exotic Z decay according to the final states and number of resonances in section IV. For each decay topology, we comment on the origin of possible UV models, provide the appropriate cuts for each topology and present the sensitivity on exotic Z decay BR. In section V, we conclude.

II. DM RELIC ABUNDANCE, INDIRECT AND DIRECT DETECTION
In this section, we briefly describe the inputs from DM direct detection, indirect detection and relic abundances employed in this study.
DM direct detection experiments look for DM collision with nuclei in the detector, which leaves visible energy in terms of phonon, electron and photon signals. We are interested in the kind of collision which provides spin-independent cross-section with nuclei, where Xenon type experiments such as XENON1T [69], LUX [70], PANDAX-II [71], provide the best sensitivity for large DM mass. For small DM mass, e.g. < 5 GeV, CRESST-II [72] and CDMSlite [73] provide better sensitivity; because their nucleus are lighter than Xenon, thus they can obtain more energy transfered from light DM collisions.
DM indirect detection experiments search for DM annihilation products like photons, electrons, positrons and anti-protons from astrophysical sources. We consider the gamma-ray line searches by Fermi-LAT [74], continuous gamma-ray limits from dwarf galaxies [75], and e ± flux measurements from AMS-02 [76]. And we consider constraints from cosmic microwave background (CMB), where DM annihilation products heat and ionize the plasma during the recombination epoch [77]. When the DM annihilation cross-section is proportional to the DM velocity square v 2 , dubbed as p-wave cross-section, the constraints from the indirect detection are negligible. If DM annihilate into two photon, the leading constraints from indirect detection normally are the gamma-ray line and CMB searches.
The relic abundance Ωh 2 = 0.12 from Planck collaboration [77] is used in this paper as a benchmark point. We assume a standard thermal freeze out. Therefore, the relic abundance only depends on the thermal average of the DM annihilation cross-section σv. 1

III. HIDDEN SECTOR MODELS AND EXOTIC Z DECAYS
In this section, we discuss several classes of well motivated dark sector models, such as Higgs portal DM, vector portal DM, inelastic DM and ALPs. These models can be probed by the exotic Z decays in future e + e − Z-factories. It is a demonstration of the capability of e + e − collider as a new physics search machine and a novel intensity frontier experiment.
For each model, we point out how it could be probed by the exotic Z decays. The existing limits from cosmology, astrophysics and collider are presented and compared with the reach of the Z-factories. If the model contains a dark matter candidate, we will derive the DM relic density by assuming thermal production. The limits from exotic Z decay are obtained from the general analysis presented in detail in section IV.4.

III.1. Higgs Portal Fermionic DM
Higgs portal is a particular simple possibility to extend the Standard Model and link it with hidden sectors. After discovering Higgs, searching for another fundamental scalar will help us improve our understanding of Electroweak Symmetry breaking. Interestingly, the other scalar is potentially related to some enigma in cosmology, such as Baryogenesis and DM. Here we will study the discovery potential of this scalar and its hidden sector by using the exotic Z decays from future e + e − collider.
The general Lagrangian of the simplified model is written down as follows [15], We assume µ 2 H < 0 and µ 2 S < 0, which trigger spontaneous symmetry breaking of the SM and hidden sector. The tree-level vacuum stability condition requires λ H > 0, λ 4 > 0; and if λ 2 < 0, |λ 2 | > λ H λ 4 /24 should be satisfied. In the broken phase, the Higgs and the singlet scalar obtain their vacuum expectation values (vevs) v H and v S , respectively, Accordingly, the DM mass m 0 χ is shifted to m χ = m 0 χ + y χ v S , which is treated as a free parameter here. Adding the extrema condition that ∂ s V = 0 and ∂ h V = 0, where V is the scalar potential, we will have the mass matrix of s and h, The scalar mass eigenstatesh ands are obtained via the following rotation, where The mass ofh ands are Let us pause here to count the relevant free parameters for the scalars. There are nine parameters including µ S , µ H , λ 1,2,3,4 , λ H and two vevs v H and v S . The extrema conditions eliminate two of them: µ S and µ H . By changing to mass eigenstate basis, the five physical observable are mh, ms, v H , v S , and mixing angle sin α, which are determined by seven parameters. Without losing generality, we set the coefficients λ 1 and λ 3 appearing in odd terms of S to be 0, which can be achieved by adding some additional quantum number or Z 2 -symmetry for S. Having observed that the Higgs mass mh = 125 GeV and v H = 246 GeV, this leads to three final free parameters ms, v S and sin α.
The decay rates and branching ratios relevant to the scalar searches are presented below. In the case that mh > 2ms, the SM Higgs decays to twos with decay width The singlet scalars can decay to pair of DM if kinematically allowed. This is the missing energy signal in the collider. The decay width is The SM Higgsh can also decay to DM pair, with a similar decay width ofs by changing cos 2 α to sin 2 α and ms to mh, In this model, the invisible decay branching ratio fors andh are The mass of the singlet scalar relevant for the study of exotic Z decays is ms m Z . If m χ < 1 2 ms, the singlet decays to DM leading to missing energy signals.
III.1.2. DM relic abundance, indirect and direct searches, and collider constraints • Relic abundance and indirect detection: In this model, the s-channel annihilationχχ →f f is the dominant process for the thermal DM freeze-out. This process is p-wave suppressed, because the mediator is CP even, while the initial state is CP odd [85]. The analytic expression for the cross-section can be written as where y f ≡ m f /v and s is the center of mass energy square. From this expression, it is clear that the annihilation cross-section is p-wave from the term s − 4m 2 χ ∝ v 2 rel . As a result, DM indirect detection can not put strong limits on this model, since the velocity dispersion of the galaxies is relatively slow. However, the temperature during DM freeze-out is relatively high, and v rel 1/3. Therefore, the p-wave suppression is not dramatic during this period.
In the DM relic abundance calculation, we consider the fermions in the final states if the annihilations are kinetically allowed. In the mass range of m b < m χ < m Z /2, the final states of quarks (b and c) and the τ lepton are included. The computation of the thermal relic density is restricted to m χ > 1.5 GeV. For the smaller DM mass, QCD non-perturbative effects and some hadronic channels should be considered. To avoid other limits, we choose m χ close to ms/2 in fig. 2. For non-resonance case, relic abundance does not lead to competitive limits.
• Direct detection: The DM χ scattering with nuclei is mediated by t-channel scalars andh, which give the possibility to detect DM via spin-independent direct detection. The spin independent scattering cross-section with nucleon is [86], where µ n is the reduced mass between DM and nucleon, f n ≈ 0.3 is the Higgs-nucleon coupling, and m n is the nucleon mass. We compare σ SI with the limits from XENON1T [69], LUX [70], PANDAX-II [71], and CRESST-II [72] as well as CDMSlite [73] for low mass DM, and show the constraints in fig. 2. The limits drop around m χ ∼ 10 GeV, because below this mass Xenon scintillators looses its sensitivity and CDMSlite becomes the dominant one.
At LEP-II, a low mass Higgs has been searched in e + e − → Z → Z * h channel, where Z decays visibly and h decays invisibly, with integrated luminosity of ∼ 114 pb −1 [61]. The Higgs bremsstrahlung process Zh is also used at higher √ s to set limit on heavier Higgs up to 114.4 GeV [94][95][96]. The searches can put constraint on sin α for the similar process Zs, which we give in fig. 2 and labeled as "LEP-Zs-inv". For the on-shell production of Zs at FCC-ee, the sensitivity on sin α has been estimated to be ∼ 0.03 for ms < 100 GeV [97]. The precision measurement of the Higgs bremsstrahlung cross-section σ(Zh) can reach the accuracy of O(0.3% − 0.7%) expected from 5 − 10 ab −1 [1,4,98], which can probe the scalar mixing down to 0.055 − 0.084 [97], labeled as "δσ(Zh)".

III.1.3. Prospects from exotic Z decay
• Exotic Z decay sensitivity: For the sensitivity at a Giga (Tera) Z-factory, we study the process Z →sZ * → (χχ) + + − , with Feynman diagram in fig. 1, wheres decays to DM particles and off-shell Z * goes to lepton pairs. We set constraints on sin α using this process and plot them in fig. 2. The previous LEP experiment [61] has searched the similar channel with Z * decay to both hadronic and leptonic channels. The details of the simulations and cuts are given in section IV.4, where the limit on the exotic decay BR has been calculated. After calculating the exotic decay BR, one can translate the Figure 1. The Feynman diagram for exotic Z decay Z →sZ * → (χχ) + + − . Note the Z is produced on shell and followed by a three-body decays + − , and the parentheses forχχ indicates they are from the decay of a resonance .   Figure 2. The 95% C.L. sensitivity for sin α from exotic Z decay Z →sZ * → (χχ) + + − at Giga (Tera) Z-factory, with y χ = 0.1(1) in the left (right) panels. We also compare with limits from DM direct detection, relic abundance, invisible Higgs BR from the LHC [87,88] (BRh inv < 0.23), the high luminosity (3 ab −1 ) LHC projection (BRh inv 0.08−0.16) [89,90] and future e + e − collider (BRh inv 0.003) [4,91] , current and future Higgs global fit from (h current global fit) [92,93] with purple and magenta lines, low mass Higgs searches in invisible channels (LEP-Zs-inv) [61,[94][95][96], and precision measurement of σ(Zh) (δσ(Zh)) [1,4,98]. The dashed (solid) lines are for existing constraints (future prospects).

LEP-Zs-inv
constraints of decay BR to physical variable sin α. We have compare our analysis with LEP and found good agreement. To be more specific, given "LEP-Zs-inv" has also worked on Z pole with an integrated luminosity 114 pb −1 , we normalize our result to the same luminosity and find the constraint is similar to the LEP. In the SM, Higgs can decay to diphoton or Zγ via top loop and W loop. Due to the mixing betweens andh, the mono-photon process Z → γs → γ(χχ) is possible. We have checked this process following the cuts in section IV.2 and found its constraint on sin α is about one order of magnitude weaker than Z →sZ * → (χχ) + + − . The main reason is mono-photon decay is loop suppressed. Furthermore, mono-photon background is higher than + − + / E background. Therefore, we do not put the constraint from mono-photon in fig. 2.
• Summary: From fig. 2, we see the relic abundance provides constraints on sin α only in the fine-tuned scenario with 2m χ ∼ ms. The indirect detection does not provide limits because it is p-wave suppressed. The direct detection provides a useful constraint, which is not sensitive to the resonant mass of ms ∼ 2m χ . At the same time, it depends on the size of the Yukawa coupling y χ . The existing and future Higgs global fit from the LHC does not provide competitive limits in comparison with precision measurement of σ(Zh), while invisible decay BR of SM Higgs provides a pretty good limit down to sin α ∼ O(10 −2 − 10 −3 ) via the existing LHC data. At the HL-LHC (3 ab −1 ), the reach of invisible BR is about 0.08 − 0.16 [89,90], which provides only a moderate improvement of the limit. The future sensitivity of BR h inv is expected to reach ∼ 0.003 at future e + e − collider [1,4,91], which can improve the limits by a factor of ∼ 8.7.
The proposed exotic Z decay Z →sZ * → (inv)+ + − can cover sin α down to ∼ 10 −2 (10 −3 ) for Giga Z (Tera Z), and such constraints do not rely much on value of y χ and χ mass. The constraints from exotic Z decay are superior than most of the existing and future searches, and only invisible SM Higgs decay search at the future Higgs-factories can provide competitive limits.

III.2. Vector portal DM
Vector portal, as another simple extension of the SM physics, employs a massive U (1) dark photon connecting the SM sector and the hidden sector [17][18][19][20][21][22]. The searches for vector portal DM and the vector field itself attract world-wide effort (see review [6][7][8] and references therein). Various experiments, such as fixed target, e + e − and pp colliders, are aiming to find such dark photon, especially utilizing its coupling to + − . Aside from decaying to SM fermions, the invisible decays of the dark photon are directly related to DM, which can be searched by radiative return process, meson decay and missing energy events in scattering processes [6][7][8].
The dark photon A , as a U (1) gauge field in the hidden sector, can mix with SM hypercharge U (1) Y field B µ through a renormalizable operator, where is the kinetic mixing parameter and c W is the cosine of the weak angle. The mass of the dark photon, m A can be obtained from Higgs mechanism in the dark sector. Interestingly, this underlying mechanism is related to our previous Higgs portal DM. We ignore here (possibly interesting) dynamics of the dark Higgs 2 . We can always rotate away the kinetic mixing terms and work in the mass eigenstate basis. The rotation is non-unitary and is written down up to where t W is the tangent of the weak angle. This formula does not apply to the region that A mass is pretty close to the mass of Z boson. In the rest of the paper, we work on the mass eigenstates of these gauge fields; without ambiguities,Ã andÃ are used to represent the mass eigenstates. After this rotation, the way that the currents couples to gauge fields are changed, and the interactions between vectors and currents up to O( 2 ) are written as follows, The massless photonÃ couples to the electromagnetic current J em . The dark photon couples to dark U (1) currents J D ; after the field rotating, a suppressed coupling to J em and J Z arises. Thẽ Z boson couples to J Z , and has the coupling to the dark currents with suppression.

III.2.1. scalar vector-portal DM
• Model: In this model, we introduce a complex scalar as DM, charged under the U (1) D , and this scalar DM interacts with the SM particles via the dark photon A . The relevant interactions can be written as follows, For m 2 S > 0 and considering Z 2 symmetry, S = 0 and S is stable. From eq. (15,16), it is clear that there is coupling between Z, A and S, which can provide interesting signal for the exotic Z decay,Z →Ã S * S from the leading terms in the Lagrangian. To have this signal, we must have this process kinematically allowed, mÃ +2m S < m Z . We will focus on the region that m S > 1 2 mÃ , such that theÃ decay dominantly to SM particles, rather than invisible DM pair. 3 The spontaneous symmetry breaking through dark Higgs φ is a simple mechanism to give mass to S and A . The difference from S is that there is no exact Z 2 symmetry to make φ stable, but the Lagrangian is similar to eq. (17), After symmetry breaking φ = 0, A and S get their mass. When φ is much heavier than S and A , it can be integrated out, and eq. (17) is enough to describe the process related to DM and various searches. When φ mass is smaller or comparable to the mass of S and A , φ need to be considered. In this case, φ can be produced at collider and decay back to 2S or 2Ã .
• DM relic abundance and indirect detection: If m S > mÃ , the dominant process controlling the freeze-out is SS * →Ã Ã . The thermal cross-section is where y ≡ mÃ /m S . By taking s → 4m 2 S , the leading term tell us that this process is s-wave. This thermal cross-section is not related to , since theÃ are produced on-shell. On the other hand, in the regime that m S < mÃ , the dominant process is S * S →Ã /Z →f f via the off-shell Z andÃ , and the thermal cross-sections for σv(SS * → ff ) are given in section VI. Since the thermal cross-section is proportional to 2 , the relic abundance will rely on the size of the kinetic mixing. This can set the target for the search of exotic Z decay. Without loss of generality, we will restrict to m S = 0.8mÃ in the parameter space, to compare various limits from the complementary experiments, shown in fig. 4. S * S →Ã /Z →f f is p-wave suppressed, which can be understood from the CP-symmetry of initial state [85]. As we discussed before, the p-wave annihilation have the suppressed signal of the indirect detection. Therefore, the corresponding limit is negligible.
• Direct detection: The scattering of S off nuclei is mediated by t-channelÃ andZ. Interestingly, the contribution fromZ exchange has been canceled by the one fromÃ coupling to J Z current [97], hence onlyÃ coupling to J em current should be considered, which can be seen directly from eq. (16). Therefore, the spin-independent scattering cross-section for S and the nucleon has a simple expression and is given below, where µ Sn = m S m n /(m S + m n ) is the reduced mass of dark matter S and nucleon n, and e is the electron charge. We add the direct detection constraints as green shade area in fig. 4.

• Existing collider limits:
Focusing on the region of mÃ < 2m S , the decay mode of the dark photon,Ã → + − , is the key channel to look for in the experiments: beam-dump, fixed target, collider, and rare meson decay. In fig. 4, we present the constraints from the experiments having the leading limits currently. There are also limits from LEP via electroweak precision observables [41]. For constraints from the LHC, the inclusive Drell-Yan process pp →Ã → + − can be used to constrain with the LHC 8 TeV data [100,101], which provides a stronger bound than the electroweak precision bounds [42,102,103]. For low mass mÃ ∼ O(GeV), the limits from B-factory is the leading one from measuring visible decay products of the dark photon, such as BaBar 2014 [104] having the limits of 10 −3 . Recently, the LHCb [105] performed dark photon search using the inclusive di-muon data. This will give the leading constraints in the mass window of (10 GeV, 50 GeV).
• Exotic Z decay search: The first process we consider is the three-body decayZ →Ã S * S → ( + − ) / E shown in the left panel of fig. 3. The limit on exotic Z decay branching ratio is given in section IV.4. Here we take the mass range ofÃ , m S < mÃ < 2m S , such thatÃ will not dominantly decay to invisible DMs, and DM relic density depends on the kinetic mixing . To constrain kinetic mixing coupling at given mÃ , we fix coupling g D , the mass ratio m S /mÃ . The corresponding limit for as a function of mÃ is given in fig. 4. The range of mÃ starts from 1 GeV. For smaller masses, other constraints like beam dump experiments become quite strong. Moreover, the exotic Z search begins to lose its efficiency due to the small separation of lepton pair fromÃ decay.
In addition to the three-body decay topology, we can also have the 2-body cascade decaỹ Z →Ã φ → ( + − )(S * S), shown in the right panel of fig. 3. This channel has resonance in both lepton pair invariant mass and invisible mass. We still consider the regime m S < mÃ < 2m S . Therefore,Ã decays into lepton pairs. We assume mφ < 2mÃ and very small Higgs portal mixing, so thatφ decay to SM particles via Higgs mixing can be neglected. Therefore, the decay branching ratio ofφ → S * S is ∼ 100%. In the right panel of fig. 4, we constrain as a function of mÃ ,   • Summary: As shown in fig. 4, LEP electroweak precision test, LHC Drell-Yan, Babar radiative return and LHCb di-muon inclusive searches can provide the direct constraints on . For mÃ < 10 GeV, Babar bounds 10 −3 , while LHC Drell-Yan and LHCb provide complementary limits 5 × 10 −3 for mÃ > 10 GeV. LEP electroweak precision test is the weakest constraint among the three.
The hint from the DM relic abundance and the constraints from direct detection and exotic Z decay rely on coupling g D . For a fixed mÃ , DM annihilation cross-section and direct detection scattering cross-section are proportional to g 2 D . The coupling for the four point vertexZ µÃ µ S * S is proportional to g 2 D , while the coupling for three point vertexZ µÃ µφ is proportional to g D mÃ . Therefore, the 3-body decay width is proportional to g 4 D , while the 2-body cascade decay width is proportional to g 2 D . For g D = 1, we see Tera Z could provide the strongest bounds at low mÃ , while direct detection provides comparable limits to exotic Z decay at high mÃ .
In comparison with the 3-body cascade decay, one might expect better constraint from 2-body cascade decay because there are resonances in both lepton pair and missing energy in this topology, while 3-body decay only has one resonance in lepton pair. This intuition is indeed correct as the sensitivity on exotic decay BR is better for 2-body cascade decay than 3-body in section IV.4, but the difference is not significant. The limits on in fig. 4 involve more parameters and couplings, which modify the dependence of mÃ .
For mÃ ∼ 1 GeV, 3-body decay loses less efficiency from lepton separation requirement than 2-body cascade decay, since the energy ofÃ in 3-body decay is generally softer than in 2-body cascade decay.
In summary, the exotic Z decay search in both topologies can provide good reach in , which is complementary and competitive to other constraints.

III.2.2. (Inelastic) vector portal fermionic DM
For vector portal fermionic DM, we consider inelastic DM model here. The constrains and future collider search of Z decays are similar to magnetic inelastic dark matter, which will be explored in section III.3.
Starting from the fermionic DM charged under dark sector U (1) D , we can write down its Dirac mass term m Dχ χ, and its Majorana mass is obtained through Yukawa interaction with a scalar Φ. The Lagrangian is The ratio of U (1) D charge of Φ and χ equals to 2. Once Φ gets vev, the DM χ gets Majorana mass along with its Dirac mass. As a result, the Dirac fermion splits itself into two Majorana fermions, which provide the DM χ 1 and its excited state χ 2 , dubbed as inelastic dark matter (IDM) [106,107].
We work with Weyl spinor and analyze the interactions for χ 1 and χ 2 . If we write χ = η, ξ † , the mass term is given as [106,108], The mass matrix can be diagonalized by a rotation, where tan 2β = 2m D /(m ξ − m η ). The mass of χ 1 and χ 2 are The vector current of the DM couples to U (1) D gauge field A . We can write both of them in the mass basis as follows, where we have define For the scalar interaction with DM can be written as, There are two interesting parameter regions for this model. In the first one, the Majorana mass is much larger than its Dirac mass, m η , m ξ m D , such that the mixing angle β is small and the mass of χ 1 and χ 2 have small corrections to its Majorana masses, where m The interactions with vector boson and scalar are mainly diagonal, while the off-diagonal interactions for χ 1 and χ 2 are suppressed.
In the second case, Dirac mass is dominant, m η , m ξ m D . Therefore, the mixing angle β is very close to its maximal value π/4. The mass of χ 1 and χ 2 are m χ 1 ≈ m D − (m ξ + m η )/2 and m χ 2 ≈ m D + (m ξ + m η )/2, with the mass splitting ∆ = m ξ + m η . We have x ≈ 2m D /(m ξ − m η ) and |x| 1. It suggests that the diagonal interactions with vector boson are suppressed while the off-diagonal interaction to χ 1 and χ 2 is dominant. For the special case of m ξ = m η , the diagonal interaction with vector boson vanishes. Since the IDM relies on the off-diagonal interactions with vector boson, the DM scattering only happens when the final states are its excited ones and provides very different phenomenology from ordinary elastic scattering in direct detection [106]. However, for scalar interactions, the diagonal interaction with the fermionic DM is proportional to m ξ + m η , while the off-diagonal interaction is proportional to m ξ − m η . The scalar mediation to diagonal terms can potentially spoil the IDM setup when the Higgs portal coupling is large.
Coming back to the exotic Z decays, we see that vector portal IDM motivates the exotic decays ofZ → χ 2 χ 1 and χ 2 χ 2 , followed by the subsequent cascade decay χ 2 →Ã χ 1 ,φχ 1 andÃ ,φ → f f , χ 1 χ 1 . It shows that IDM with vector portal can motivate the topologies of exotic Z decay in section IV.

III.3. Magnetic inelastic DM, Rayleigh DM
The coupling of DM to the Standard Model particles can be very weak. One possible scenario is that the hidden sector interacts with the Standard Model via high dimensional operators. The representative models, the Magnetic inelastic DM (MIDM) and Rayleigh DM model (RayDM) [36][37][38][39][40], are introduced, and their relevance to the exotic Z decay is studied in this section.
The two models, the MIDM and RayDM, can be derived from the same UV model [40], χ is fermionic DM with a Dirac mass term m χ and Majorana mass term δm. It interacts with scalar φ and another fermion ψ via a Yukawa coupling. The Dirac and Majorana mass terms can split DM χ into two Majorana fermion χ 1 and χ 2 , where we assume m χ 2 > m χ 1 . The fermion ψ and scalar φ have the same charge under SM gauge group SU(2) L × U(1) Y [40]. The dark matter will couple to photon via ψ and φ loop. Integrating out ψ and φ will generate two higher dimensional operators. The first operator is MIDM operator [36][37][38], and the second is RayDM operator [39]. Both of them are given below, Note there are also operators including γ 5 in the DM bilinear, which corresponds the electric dipole operator. For RayDM, the corresponding one is whereB µν = µναβ B αβ and µναβ is the anti-symmetric Levi-Civita symbol. The interaction scale Λ has been calculated in [40] 1 where we have assumed that ψ and φ are singlet under SU (2) L and are charged under U (1) Y . In eq. (31), we have assumed φ mass is similar to M ψ and we take the form factor function to be O(1). These two operators can lead to the cascade decay Z → χ 2 χ 1 → (χ 1 γ)χ 1 and the three-body decay Z → χ 1 χ 1 γ at Z-factory, with Feynman diagrams given in fig. 5. In the exotic Z decay study, we will choose a significant mass splitting between χ 1 and χ 2 , to get a hard photon signal which can be detected at Z-factories.

III.3.2. DM relic abundance, indirect and direct searches, and collider constraints
• Relic abundance and Indirect detection: We focus on the case that there is a significant mass splitting between χ 2 and χ 1 , which can give rise to interesting photon signal in exotic Z-decay. In this case, the relevant annihilation initial state contains only χ 1 . The annihilation rate is dominated by the Rayleigh operator into γγ, γZ, ZZ and W + W − . For the mass range m χ 1 < m Z , we find only the following annihilation cross-section relevant [39], where y ≡ m χ 2 /m χ 1 . The two annihilation cross-section for RayDM are for O RayDM and O γ 5 RayDM respectively, where the former one is p-wave suppressed while the second one is s-wave. The annihilation process for the MIDM scenario is two loop suppressed. This can be seen in eq. (32), eq. (33) and eq. (34) through the dependence on Λ MIDM and Λ RayDM , respectively. The annihilation into gamma-ray lines is constrained by Fermi-LAT search [74] (blue shaded region ) and also by CMB [77] (purple shaded region), which we constrain M ψ as a function of m χ in fig. 6. The long dashed lines are for O γ 5 RayDM , while the dashed lines are for O RayDM which is very weak due to the p-wave suppression.
• Direct detection: In the case of large splitting, only Majorana χ 1 is relevant for direct detection, because inelastic scattering into χ 2 is kinetically forbidden. Therefore, the scattering cross-section is dominated by the loop exchange of two photons from Rayleigh operator, and the spin-independent cross-section per nucleon is given below [39], where m N is the mass of nuclei N , A is the nucleon number, Z is the proton number of nuclei and Q 0 is the nuclear coherence scale Q 0 = √ 6(0.3 + 0.89A 1/3 ) −1 fm −1 . The current leading constraints on spin-independent cross-sections are XENON1T [69], LUX [70], PANDAX-II [71], and CRESST-II [72] as well as CDMSlite [73]. The limits from direct detection constraints are shown as magenta in fig. 6. Only O RayDM is shown as in dashed lines, because O γ 5 RayDM produces spin-dependent cross-section and spin-independent cross-section is suppressed.
• Existing collider constraints: Besides DM indirect and direct detection, the MIDM and RayDM operators can also get constraints from mono-jet and mono-photon searches at LHC and LEP.
The Rayleigh operator O RayDM has been studied in mono-photon, mono-jet and mono-V searches [110], where V stands for vector gauge boson W and Z. The authors found that the limits from mono-photon provides the strongest bound and constrain Λ RayDM 510 GeV at 95% C.L, for m χ 1 100 GeV from the LHC 8 TeV at 20 fb −1 [111,112]. Very recently, ATLAS [113] has explored 13 TeV data to search mono-photon signature with integrated luminosity 36fb −1 , and it pushes the limit to Λ RayDM 725 GeV. These limits has been integrated in the right panel of fig. 6, and denoted as "mono-γ". For O γ 5 RayDM , the limits are similar as O RayDM and therefore we only show the results for O RayDM .
For the MIDM operator O MIDM , [109] has studied the limits from the mono-jet, mono-photon and di-photon searches at 8 TeV and 14 TeV LHC. For a significant splitting, they found monophoton search [111] is the most stringent, similar to the RayDM operator case. For m χ 1 = 10 GeV, it requires Λ MIDM 2400 GeV, and the result is roughly unchanged for m χ 2 > 20 GeV. In the left and middle panels of fig. 6, we vary DM mass m χ 1 from 0 to 40 GeV. Since its mass is much smaller than the required photon p T and MET, we expect the constraint to be similar as m χ 1 = 10 GeV. For mono-photon search at the LHC 14 TeV with 300 fb −1 , the corresponding limit is estimated to be Λ MIDM 8200 GeV [109], and labeled as "mono-γ" in fig. 6. Figure 5. The Feynman diagrams for the cascade decay process Z → χ 2 χ 1 → χ 1 χ 1 γ from O MIDM and the three-body process Z → χ 1 χ 1 γ from O RayDM .
Fe rm i-L AT CM B  Figure 6. The 95% C.L. sensitivity for M ψ from exotic Z decay Z → / E + γ, for MIDM operator in the left (middle) panels with different mass splitting and for Rayleigh operator in the right panel. The constraints are labeled as Giga Z and Tera Z for future Z-factory with λ = 4π, and the LEP limit from [64] is shown. We also compare the limits from DM direct detection, indirect detection searches, mono-photon, and mono-jet searches at the LHC. For RayDM, the gamma-ray constraints from Fermi-LAT and CMB use long dashed line for O γ5 RayDM and dashed line for O RayDM . For collider limits, the two operators are similar and for spin-independent direct detection limits, only O RayDM is constrained.
For the MIDM case, it is interesting to note that, when m χ 2 = m χ 1 , the exotic Z decay Z → / Eγ loses its sensitivity at Z -factory, and also for mono-photon search at the LHC. The mono-jet search will be better than the mono-photon search in this case. Moreover, [109] pointed out that actually the invisible decay width measurement of Z can beat the mono-jet search at the LHC 14 TeV with 3 ab −1 integrated luminosity, which suggest M ψ 226 GeV for m χ 1,2 = 10 GeV. We have plotted the invisible Z width constraint in panel (a) of fig. 6.
Given the high center of mass energy at the LHC, it can search for the EW charged particles ψ and φ directly from Drell-Yan production and their subsequent cascade decays [114]. The Drell-Yan search could be more restrictive than mono object searches, but this conclusion is very model dependent, see [114]. For example, when ψ and φ are SU (2) L singlet, or they decay dominantly to tau lepton and (or) gauge bosons, the sensitivity from Drell-Yan is very poor, even at the LHC 14 TeV with 300 fb −1 .
For mono-photon at LEP, the L3 collaboration has collected data with 137 pb −1 at the Z pole, which can limit the BR of exotic decay Z → γ / E down to 1.1×10 −6 if photon energy is greater than ∼ 30 GeV [64]. The OPAL collaboration has a similar study at Z pole but with only 40.5 pb −1 [65]. There are also many off-Z peak measurements on single photon final state. The one with 176 pb −1 data taken at 189 GeV has been carried out by the L3 collaboration, which looks for MIDM topology Z → χ 2 χ 1 → (χ 1 γ)χ 1 , and bounds the cross-section of such topology to be smaller than 0.15 − 0.4 pb with some dependence on m χ 1 and m χ 2 [66]. The leading constraint is from L3 measurement at Z pole due to large resonant cross-section, and we label the constraints as "L3 Mono-γ" in fig. 6. We see that this constraint is comparable to the future LHC reach in middle panel of fig. 6.

III.3.3. Prospects from exotic Z decay
• Exotic Z decay sensitivity: For exotic Z decay with final state / Eγ, we summarize the results of the cascade decay process Z → χ 2 χ 1 → χ 1 χ 1 γ from O MIDM and of the three-body process Z → χ 1 χ 1 γ from O RayDM in fig. 5. The limits on such exotic decay BR is given in section IV.2, and we can calculate the limits for Λ MIDM and Λ RayDM accordingly, then convert them into constraints for M ψ by eq. (31). The limits are given in fig. 6, and labeled as "Giga Z" and "Tera Z".
The results of the MIDM operator are presented in the left and middle panel of fig. 6. We find exotic Z decay can reach M ψ ∼ O(10 4 ) GeV, which is much better than mono-photon searches at the HL-LHC with M ψ ∼ 10 3 GeV. The production cross-sections for χ 2 χ 1 at Z-factory and LHC both scales as 1/Λ 2 MIDM . However, the cross-section at Z-factory benefits from Z resonance comparing with at the LHC, therefore have larger statistics. Moreover, the / E + γ searches at Zfactories have much cleaner environment than hadron collider. As a result, exotic Z decay can give the M ψ reach two orders better than mono-photon search at the LHC or HL-LHC. The indirect detection of gamma lines at Fermi-LAT provides a similar constraint to the LHC 8 TeV. The direct detection does not provide any constraint for MIDM operator, because the mass splitting between χ 1 and χ 2 is too large.
In the right panel of fig. 6, for RayDM operator, we find the mono-photon search at the LHC can reach M ψ to a few hundreds of GeV, which is better than exotic Z decay with M ψ 100 GeV. The reason is that the cross-section for χ 1 χ 1 γ is proportional to s 2 /Λ 6 RayDM ,. Since the Z-factory has a small center of mass energy square s ∼ m 2 Z , it has less sensitivity. The constraint from direct detection is very weak, because it is a two loop process. The gamma line constraint from Fermi-LAT is comparable to other constraints and is strongest at m χ 1 around 100 GeV.
• Summary: We find complementarity between exotic Z decay Z → / Eγ at Z-factory and mono-jet or monophoton search at the LHC with large mass splitting between χ 1 and χ 2 . For very small mass splitting, the photon from cascade decay χ 2 → χ 1 γ becomes very soft, thus the mono-photon and mono-jet search via initial state radiation are better. However, invisible Z width measurement can provide a better limit M ψ 226 GeV. For MIDM operator, future Z-factory can provide the leading constraints, while for RayDM operator, the HL-LHC can provide better constraints.
III.4. Axion-like particle Axion-like particle (ALP) is a light pseudo-scalar which couples to gauge fields via anomalous terms and interacts with fermions with derivatives, ∂ µ aψγ µ ψ. Its presence is quite generic in UV theories, such as string theory [30,32,34], and Supersymmetry [26][27][28]. It can be a portal connecting dark matter with the standard model sector [31], and ultralight ALP is dark matter candidate by coherent oscillating in the universe [115][116][117]. Recently the dynamics of ALP in the universe has also been proposed to solve the Higgs hierarchy problem [118]. For our Z-factory study, we are focusing on the mass range of ALP from 0.1 GeV to Z boson mass. Although we focus on the case of ALP, our analysis and results in this section can be applied to scalar easily.   Figure 8. The limit on Λ aBB , ALP coupling to hypercharge field, from future Z-factory. The limits from LEP I [119] γγ search, LEP II (OPAL) 2γ and 3γ searches [120], , LEP (L3) 3γ searche at Z pole [67], ATLAS 3γ and Z → 3γ [121,122] search are translated to limits on Λ aBB following [123]. There are three type of signals Z → 2γ, 3γ and / Eγ, depending on m a . In / Eγ final state where a decay outside the detector, we have set the detector length to be 6 meter and LEP limits on this final state from L3 collaboration [64] has been plotted.
ALPs can have interactions with standard model particles fermions, gauge fields, Higgs obeying the (discrete-)shift symmetry. Here, we focus on the ALP coupling to the U (1) Y gauge field B µ 4 , This interaction gives the decay rate of the ALP as Γ(a → γγ) = 1 64π and the rate of the Z decay, Depending on the a → γγ decay length, the analyses are performed in the two separate regimes: one is ALP decaying inside the detector, and the other is decaying outside the detector. For decay inside the detector, we focus on the prompt search, and leave the interesting case of displaced vertex to future work. For decay outside the detector, the signal is mono-photon + / E. The transverse radius of the detector radius is taken to be 6 meters. The decay length of the ALP is computed according to the boost γ a of the ALP, D ≡ γ a cτ a , where the γ a = E a /m a is the boost and τ a = 1/Γ a is the lifetime of a. Since the initial state is Z boson at rest and the final state is aγ, the energy E a is fixed by m a . D = 6 m is plotted in fig. 8 as a dotted black line. Below it, the ALP has a decay length D smaller than 6 m. However, it can still decay outside the detector with a probability of 1 − e −D/(6 m) . We account for this probability to rescale the signal events in the detector, which leads to sensitivity below the line. In the prompt decay region, for the high mass axion, the boost of axion is small, the dominant channel to search for ALPs is 3γ. When the mass of the ALP is below O(1) GeV, the boost of axion makes the two photons from axion decay close to enough, and cannot be resolved. The 2γ search channel is more relevant.
The current constraints for this operator are given by LEP and LHC photon searches. In fig. 8, the LEP I [119] uses inclusive di-photon search e + e − → 2γ + X covering the small mass region. In the higher mass region, the boost of the axion decreases and 3γ channel is considered. The LEP II (OPAL) have 2γ and 3γ data [120], which are employed to put the bounds on the process, e + e − → γ/Z → aγ → 2γ + γ. The L3 collaboration has searched the process Z → aγ → (γγ)γ at Z pole, with limit on BR of order 10 −5 [67]. ATLAS 3γ and Z → 3γ [121,122] search can be translated to the ALP bound as derived in [123].
For / E + γ search, the strongest bound from LEP comes from L3 collaboration with 137 pb −1 data at the Z pole [64] as discussed in section III.3. It can limit the BR of exotic decay Z → γ / E down to 1.1×10 −6 if photon energy is greater than ∼ 30 GeV. It directly excludes Λ aBB < 4.3×10 4 for Z → / E + γ decay, and we label it as "L3 ( / Eγ)" in fig. 8. In the Z-decay search, the ALP will give topologies Z → / E + γ and Z → 3γ, 2γ, depending on the life-time and boost of the ALP. Z-factory limits on the ALP are given in fig. 8, which is about two order of magnitude better than the current constraints from LEP and LHC.

IV. SEARCHING FOR EXOTIC Z DECAYS AT FUTURE Z-FACTORIES
In this section, we make projections for the sensitivity of exotic Z decay searches at future Z-factories. Motived by the previous discussed dark sector models, we classify decay channels by final states, the number of intermediate resonances, and different topologies. In most of the cases, we clarify the connections between the potential models and each topology. As Z is neutral, the final states of its decay can be described as Since lepton and quark are charged, they will show up in pairs. The n is referred to as the number of particle or pair of particles. In our analysis, we choose to consider the number of final state particles to be less than 5. The / E can be considered as two particles, since normally it is constituted of two DM particles. It also can be a neutral particle which does not interact with detector and decays outside of it. The final states can be further grouped according to whether they are the decay products of some intermediate resonance. This resonance can be the mother particles for (γγ), + − , (qq) and / E. The kinematic information of the resonance decay can help us improve the search strategies. The details of classification are given in table I. The first set of channels has the missing energy in the final states. Since electron collider has full kinematic information of initial states, the missing 4-momentum can be fully reconstructed. This is the major advantage of electron collider compared with hadron collider in searching for exotic Z decay with missing energy. The second set of channels does not include missing energy. They are pure jet final states (jj)(jj), (jj)(bb), (bb)(bb) and three photon final state γγγ. They can come from dark sector particles decays, which do not involve dark matter. Due to the cleaner environment of electron collider, it is better than hadron colliders to measure pure hadronic final states. For jjjj final state, since it has large SM background, we concentrate on the case where it has two resonances. When generating corresponding SM backgrounds, one additional photon is included to count the initial state radiation (ISR). The on-shell intermediate particles should be neutral, since LEP searches have already put severe constraints on charged particles with mass smaller than m Z /2.
In the following subsections, we will discuss the possible models and the sensitivity of each channel at future Z-factory. The section IV.1 introduces the basic setup and performance for future Z-factories at FCC-ee and CEPC, and explores the sensitivity of exotic Z BR at this future Z-factory for different topologies from section IV.2 to section IV.7. To compare the future Z-factory and HL-LHC, section IV.8 presents the reach on those exotic Z BR for the HL-LHC. The summary of this comparison between the future Z-factory and HL-LHC is in fig. 16.

IV.1. Performance of Future Z-factories
The exotic Z decay phenomenology at future Z-factories at studied in this section. A Z-pole run has been considered for both FCC-ee and CEPC [124,125]. Given that the measured cross-section of hadronic Z is 30.5 nb [5], the integrated luminosity for Giga Z (10 9 Z) and Tera Z (10 12 Z in the plan of FCC-ee) are 22.9 fb −1 and 22.9 ab −1 , respectively.
We simulate the backgrounds and signals in the electron-positron colliders at the Z mass energy using MadGraph5 aMC@NLO [126] and analyze them at parton level. Assuming that the detector performance is similar for different future electron colliders, we follow the detector effects at CEPC [4] and apply the following Gaussian smearing in our analysis: Photon energy resolution: exotic decays topologies n res models    Lepton momentum resolution: ∆ GeV Jet energy resolution: We make conservative assumptions about the tagging efficiency: 80% for about b-tagging efficiency, 9% for c quark mis-tagging rate and 1% for light flavor mis-tagging rate [4]. We also require that all visible particles satisfy |η| < 2.3 (cos θ < 0.98). In addition, the photon, lepton and jet energy should be larger than 10 GeV. For events with missing energy, we require / E > 10 GeV as well. Lastly, both the photons and electrons in final state are separated by θ ij 10 • = 0.175 radian. The charged leptons normally have better resolution than photons, thus the separation requirement that we choose here is conservative. For jets, we use a conservative separation requirement θ ij 0.4 radian corresponding to ∆R ≥ 0.4 at LHC 5 . The study for LEP3 (a 240 GeV circular ee collider using LHC tunnel) with the CMS detector [128] shows the jet angular resolution can be 30 milliradian for energies below 100 GeV. The separation requirement for jets at lepton collider could be optimized due to much less QCD backgrounds than LHC in principal. We leave the optimization for lepton collider as the future study.
To derive the exclusion limits, the confidential level for the sensitivity calculation adopts Poisson probability [129]. When background event number B 1 , the significance is about S/ √ B which is proportional to √ L, where L is the integrated luminosity. Therefore, the sensitivity reach of Giga Z and Tera Z differ by about 10 1.5 . When background event number B 1, the Poisson distribution with zero background assumption leads to a constant limit for signal. In this case, the exclusion limit is linear to L, thus Giga Z and Tera Z differ by about 10 3 . If B 1 for Giga Z while B > 1 for Tera Z, the difference of the sensitivity reach is in the range of 10 1.5 − 10 3 .
Te ra Z BR(Z→γ+a) Figure 9. The 95% C.L. exclusion on exotic Z decay BR for the final state Z → / Eγ. (a): the decay topology 1A, Z → χ 2 + χ 1 → χ 1 γ + χ 1 from MIDM model. The numbers in each block are the sensitivity reach for the exotic Z decay BR in log 10 for Giga Z and Tera Z respectively, while the color mapping is coded for Tera Z. (b): the decay topology 1B, Z → χχγ from RayDM model. (c): the decay topology 1C, Z → aγ → / Eγ from axion-like particle model.
In this section, we discuss the exotic Z decay with the final state / E + γ. We consider the decay topologies Z → χ 2 χ 1 → χ 1 γ +χ 1 and Z → χχγ, where χ and χ 1 are fermionic DM. χ 2 is an excited DM state which decays back to χ 1 . We also consider 2-body decay Z → aγ → ( / E)γ, where a is a pseudo-scalar as missing energy signal if it is stable at collider scale or it decays to dark matter particles. We denote these three topologies as 1A, 1B and 1C, respectively, shown in table I. The UV models for the 1A and 1B are the MIDM and RayDM model, while 1C is motivated by ALP. The fourth topology, denoted as 1D, is Z → A γ → (χχ)γ. It can come from the Wess-Zumino term µνρσ A µ B ν ∂ ρ B σ , when the dark photons couple to anomalous currents [130,131]. After integrating by parts, the longitudinal part of A have similar interaction as the topology 1C; thus the limit on exotic Z decay BR is similar as 1C.
The SM backgrounds for these final states are mainly e + e − → γνν. In our simulation, we include one more photon to account for the ISR effect. For γν eνe , this process is mediated by either off-shell W ± * or off-shell Z * , while γν µνµ and γν τντ are mediated by off-shell Z * . In these processes, most of the γs come from ISR, or internal bremsstrahlung via the t-channel W boson. The background photons are generally quite soft due to their origin as ISR.
The three models have the different kinematic distributions for the mono-photon. For the topology 1A, Z → χ 2 + χ 1 → χ 1 γ + χ 1 , the photon energy spectrum has a box shape due to the cascade decay in this process. The minimum and maximum of the photon energy are The distribution of photon energy is flat between E min,1A γ , E max,1A γ , and the edge of photon energy distribution can be used to determine the mass of DM. Therefore, aside from the pre-selection cuts, we further impose the cuts below, where m inv is the invariant mass of missing energy. The second cut comes from momentum conservation that the invariant mass of a set of particles is larger or equal to the sum of individual masses. According to the recoil mass relation, E γ and m inv are not independent with each other. If we apply the first cut, the second cut is automatically satisfied. Nevertheless, we list the second cut, sine this is not redundant in other cases.
For the topology 1B, Z → χχγ, it has a broad distribution in photon energy. The recoil mass m inv is related to the photon energy E γ by In the mean time, the relation m inv ≥ 2m χ gives the maximum allowed photon energy Thus, in addition to the pre-selection cut, we imposes the following cuts to further suppress the SM background, 1B : The lower bound of E γ is chosen to keep significant amount of signal event, and to reject SM background as much as possible.
For the topology 1C, Z → aγ → ( / E)γ, the photon energy spectrum is a delta function with Considering the photon energy ∼ 10 GeV, the energy resolution for this photon energy is around 5% according to eq. (40). Therefore, we can choose a 2 GeV window on the photon energy, After applying the pre-selection cuts and the specific cuts for the topologies 1A, 1B and 1C, we obtain the 95% C.L. exclusion on the exotic Z decay BR in fig. 9. In the panel (a) of fig. 9 for the topology 1A, the numbers in each block are log 10 (BR) for Tera Z (black) and Giga Z (dark red). It is clear that the sensitivity on BR for Giga Z and Tera Z differ by a factor of 10 1.5 . The reason is SM background γνν has event number much larger than 1 for both Tera Z and Giga Z, thus the sensitivity is scaled as S/ √ B. As a result, the sensitivity scales with luminosity as √ L, so the BR sensitivity gets a factor of 10 1.5 increase from Giga Z to Tera Z. For Giga Z, the limit on BR falls in the range 10 −6 − 10 −7 , while reaches 10 −7 − 10 −8 for Tera Z. In the panel (b) of fig. 9 for the topology 1B, the luminosity scaling between Giga Z and Tera Z is the same as in 1A. The limits on BR for Giga Z is close to ∼ 10 −6 , which is a little bit weaker than 1A due to its 3-body decay topology. In the panel (c) of fig. 9 for the topology 1C, the luminosity scaling between Giga Z and Tera Z is similar to 1A and 1B. However, the sensitivity on BR for Giga Z is close to ∼ 10 −8 , which is about 2 orders better than 1A and 1B. The massive resonance in / E implies that the energy of photon is mono-chromatic, which greatly reduces SM background.
In this section, we focus on the exotic Z decay to the final states / E + γγ. The decay topologies can be classified by number of resonances. SM background for this final state is coming from e + e − → γγνν. The general feature of the background is the same as γνν, where the photons dominantly come from ISR and tend to be soft.
For topologies with 2 resonances, the first one is the topology 2A, Z → φ d A → (γγ)(χχ), where A is a vector boson which decays into a pair of DM and φ d is a scalar which decays into a pair of photons. It can be motivated by the vector portal model in section III.2.1. The dark Higgs φ d decays to diphoton via SM Higgs mixing, or by the loop of heavy vector-like charged particles. The dark photon decays to fermionic DM which is charged under this U (1) .
The second topology with 2 resonances is the topology 2B. Z → φ A φ H → (χχ)(γγ), where φ A and φ H are CP odd and CP even scalar respectively. The topology 2B can be well motivated by the two Higgs doublet model (2HDM). The CP even scalar φ H is the mixture of CP even scalars in the 2HDM, and can decay to diphoton via loop. For the CP odd scalar φ A , decaying toχχ, one needs to add a singlet CP odd scalar φ a which couples to DM via iφ aχ γ 5 χ. The φ a can further couples to scalars by iφ a H † 1 H 2 + h.c. [132], where H 1,2 are the doublet Higgs in 2HDM. After working out the mass eigenstate, φ A is the mixture of singlet CP odd scalar φ a and doublet CP odd scalar in H 1,2 . As a result, it can have the decay topology Z → φ A φ H and φ A can further decay toχχ. Since the topology 2A and 2B has the same kinetic feature, the sensitivity to them are similar. Due to the similarity, we take the topology 2A as an example. With the presence of two resonances in γγ andχχ, we propose to use the following cuts besides the fiducial selection, Note our invariant mass window cut for diphoton γγ and missing energyχχ are conservative. The resolution for diphoton invariant mass is about 0.5 GeV at LEP [120]. The invariant mass of missing energy is determined by the energy resolution of the diphoton system, which should be smaller than 2 GeV according to eq. (40).
For the topology with 1 resonance, we have the topology 2C, Z → χ 2 χ 1 , with the subsequent decays of χ 2 → χ 1 φ d → χ 1 (γγ), where χ 1,2 are the light and heavy DM, and φ d is a scalar. This topology can be realized by either MIDM model in section III.3 or IDM embedded in vector model in section III.2.2. Since there is a resonance in γγ, one can propose the following cuts besides the pre-selection cuts, For the topology with 0 resonance, we have 2D, Z → χ 2 χ 2 , with the subsequent decay of χ 2 → χ 1 γ. This topology can be motivated by an extended MIDM model as explained in section III.3.1. From the event topology, the two photons in final state has no resonance feature. However, the photon energy distribution has a box shape similar to model 1A. The topology dictates the energy range of both photons, Therefore, we propose the following cuts besides the pre-selection cuts for model 2D, In fig. 10, we show the 95% C.L. exclusion on exotic Z decay BR for the final state Z → / Eγγ. In panel (a) of fig. 10 with two resonances and m χ 1 = 0 GeV, the SM background events are so suppressed that the event number is typically smaller than 1. As a result, the sensitivity does not scale as S/ √ B ∼ √ L, that the scaling is linear with ∼ L. This behavior can be seen from panel (a) in fig. 10 that BR sensitivity of Giga Z is around [10 −8.4 , 10 −6.7 ] and [10 −11 , 10 −9.7 ] for Tera Z. The sensitivity difference of these two are about 10 −2 − 10 −3 . The best sensitivity appears near the region where m A ≥ 10 GeV because of the pre-selection cut E inv > 10 GeV. The sensitivity gets better when m φ d becomes large since the photon becomes more energetic and the SM background becomes smaller. In panel (b) of fig. 10, we assume m χ 1 = 0. With only one resonance, one should expect the sensitivity of figure (b) to be weaker than the sensitivity in figure (a) with two resonances. We do see this point that sensitivity for figure (a) is better than figure (b) at the same scalar mass m φ d = m φ for Tera Z case, but not for Giga Z. We have looked into the cut efficiency of signal and background, which explains such behavior. The cut efficiencies for signal in respectively, which shows that the resonance condition for missing energy does help to reduce the SM background. In Giga Z case, the background event is already smaller than 1 for figure (b), therefore it has slightly better sensitivity than figure (a) due to higher signal acceptance. For the case of Tera Z, the increased luminosity has brought back the needs to suppress the SM background, therefore figure (a) has better sensitivity than (b). In panel (c) of fig. 10, the limits on exotic Z decay BR is not as good as figure (a) and (b), because there is no resonance feature in the topology. However, the constraints can still reach [10 −8.4 , 10 −7.4 ] for Giga Z and [10 −10.3 , 10 −9.2 ] for Tera Z.
For the panels (a), (b) and (c) in fig. 10, one might expect the sensitivity on BR decreases because the number of resonance n res decreases. This is clearly true when comparing n res = 1, 2 with n res = 0. However, for n res = 2 and n res = 1, the difference in sensitivity is not very significant, while the sensitivity relies more on the particle mass and the cuts. For example, the best sensitivity for 2A appears when m A ∼ 15 GeV and m φ d ∼ 60 GeV. The higher m φ d the higher photon energy, however one should also keep m A large enough to pass the missing energy cut / E > 10 GeV. The best sensitivity for 2C appears when m χ 2 ∼ 90 GeV and m φ d ∼ 80 GeV if fixing m χ 1 = 0 GeV. This high m φ d mass can guarantee a harder photon spectrum than 2A. Therefore, even without the resonance cut on / E, the SM background of 2C is similar to that of 2A, making the sensitivities on BR are similar.
In this section, we focus on the exotic Z decay to final state / E + + − . The SM background for this final state is coming from + −ν ν, mediated by off-shell gauge boson γ * , Z * and W * . Comparing with ISR photon, the energy spectrum of leptons are harder. And the spectrum of invariant mass of + − is softer than that ofνν. Given the fact that when mediated by W * , and ν are sharing similar kinetic distribution, this can not lead to difference in invariant mass. The difference is originated from that + − can be produced from γ * favoring smaller invariant masses, whileνν from Z * has much harder spectrum because small m inv are suppressed by a factor m 4 inv /m 4 Z .  Figure 11. The 95% C.L. exclusion on exotic Z decay BR for the final state Z → / E + − . The numbers in each block are the sensitivity reach for the exotic Z decay BR in log 10 for Giga Z and Tera Z respectively, while the color mapping is coded for Tera Z. The dark Higgs bremsstrahlung process can be naturally realized by vector portal model with a dark Higgs in section III.2.1.
The topologies with 1 resonance are 3B, 3C and 3D. The topology 3B is a 3-body decay Z → A S * S → ( + − ) / E, which can be motivated from vector portal model with scalar DM in section III.2.1. The topology 3C is also a 3-body process mediated by off-shell Z or photon, Z → φ(Z * /γ * ) → ( / E) + − where φ is assumed to decay outside of detector. It can be motivated by axion-like particle model in section III.4, or Higgs portal model in section III.1 where φ is a singlet scalar which mixes with SM Higgs and can decay to DM pairχχ. The topology 3D is a 2-body cascade decay, Z → χ 2 χ 1 → A χ 1 + χ 1 → ( + − ) + / E, which can be motivated by Vector portal and Inelastic DM in section III.2.2.
We will study the constraints from exotic Z decay in topologies 3A, 3B, 3E and 3F. They are chosen to represent different n res and number of particles in the cascade decay, from 2-body to 4-body. Besides the same pre-selection cuts, we propose the following different cuts for different topologies: In fig. 11, we show the constraints on exotic Z decay branching ratio BR(Z → / E + − ). For Giga Z, the topologies with n res > 0 will probe exotic Z decay BR down to ∼ 10 −8.5 , while for n res = 0 the sensitivity of BR can reach ∼ 10 −8 . With more resonances, the SM background events are suppressed a lot that the event number is typically smaller than 1. As a result, the sensitivity reach scales as L. For n res > 0 in panels (a) and (b) of fig. 11, the sensitivity on BR between Giga Z and Tera Z differs by factor of 10 2 → 10 3 due to small number of SM background. While for n res = 0 in panels (c) and (d), the sensitivity on BR between Giga Z and Tera Z differs by factor of ∼ 10 1.5 which is a very typical scaling from S/ √ B ∼ √ L.

IV.5. Z → / E + JJ
In this section, we focus on the exotic Z decay to final state / E + JJ. The J includes both the light flavor jets j and bottom quark jets b. The topologies are similar as Z → / E + + − , and the limits on exotic Z decay BR are calculated through the same procedure. The SM background for this final state is dominantly fromνν + JJ, mediated by off-shell gauge boson γ * , Z * and W * .
We choose three topologies 4A, 4B and 4C to study the sensitivity reach of exotic Z decay BR. The topology 4A is Z → φ d A → (χχ)(jj), and 4B is Z → φ d A → (bb)(χχ). Both topologies can be motivated by the vector portal model in section III.2. Here we do not use φ d → jj because Yukawa coupling is suppressed light quark mass. The last topology 4C is Z → χ 2 χ 1 → bbχ 1 + χ 1 → bb / E, which can be motivated from the MIDM operator in section III.3.1. Besides the fiducial cuts, we propose the following cuts for different topologies: In fig. 12, we show the constraints on exotic Z decay branching ratio BR(Z → / EJJ). For Giga Z, the exotic Z decay BR can be probed down to 10 −7 − 10 −8 , while the sensitivity of Tera Z is generally better by factor of ∼ 10 1.5 comparing with Giga Z. Comparing the BR sensitivity of 4A and 4B, we see that the difference between light flavor jet j and heavy flavor jet b is not large. One might expect the sensitivity of / E(bb) should be better than / E(jj), due to smaller SM background. However, the topologies 4A and 4B are not the same, where in 4A the jets come from A while in 4B the b-jets come from φ d .
In this section, we focus on the exotic Z decay to final state (JJ) + (JJ). Note that we only discuss the cases where there are two jet resonances in the final states. The SM background for this final state are mostly from electroweak process, mediated by off-shell gauge boson γ * , Z * and W * .
In table I, we have listed the topologies. 5A could be motivated from Higgs bremsstrahlung in vector and scalar portal model. We will choose the topology 5A to illustrate the sensitivity to the BR of (JJ) + (JJ) final state. We divide the final states with three combinations (jj) + (jj), (jj) + (bb) and (bb) + (bb), where the last two are denoted as 5B and 5C. There could be other topologies like Z → φ A φ H → (JJ)(JJ) from 2HDM, but the topology and kinematics are similar, therefore their sensitivity should be similar to 5A. Beside the pre-selection cuts, we add the following similar cuts for the topology 5A, 5B and 5C: Figure 13. The 95% C.L. sensitivity for exotic Z decay Z → (JJ)(JJ), where J could be light flavor jet or b-jet. The numbers in each block are the sensitivity reach for the exotic Z decay BR in log 10 for Giga Z and Tera Z respectively, while the color mapping is coded for Tera Z. The decay process is Z → φ d A with subsequent decays φ d → jj and A → jj. We show three combination (jj)(jj), (jj)(bb) and (bb)(bb) in the figure.
5B : The χ 2 method are employed to determine which pair of jets are from A decay or φ d decay. The mass window that we take is conservative. For example, at E j = 40 GeV, the jet energy resolution is about 5% leading to ∆E j = 2 GeV from eq. (42). In fig. 13, we show the constraints on exotic Z decay branching ratio BR(Z → (JJ)(JJ)). For Giga Z, the exotic Z decay BR can be probed down to ∼ 10 −5 for (jj)(jj) final state, ∼ 10 −6 for (jj)(bb) and 10 −6.5 for (bb)(bb). The sensitivity of Tera Z is generally better by factor of ∼ 10 1.5 comparing with Giga Z, from the integrated luminosity scaling S/ √ B ≈ √ L. It is clear that the sensitivity for heavy flavor jet is slightly better than light flavor jet. This is because the heavy flavor jet has fewer SM background events, by a factor of N 1/2 f ≈ 10 0.5 , where N f is the number of flavor in jets.

IV.7. Z → γγγ
In this section, we discuss the exotic Z decay to final state (γγ)γ. The SM background for this final state γγγ are dominated by QED process e + e − → γγ with an extra γ from initial state radiation, therefore the photon energy generally tends to be soft. The signal topology 6A in table I, Z → φγ → (γγ)γ, could be motivated from axion-like particle, or from Higgs portal scalar, which can decay to γγ from top loop. We take axion-like particle as an example in fig. 14, and the result should also apply to Higgs portal scalar. Besides the pre-selection cuts, we propose the following cuts for topology 6A where E 6A γ = (s − m 2 φ )/(2 √ s). We use the χ 2 method to determine the pair of photons from φ decay and single out the 3rd photon. The energy of the 3rd photon E 3rd γ is very close to E 6A γ , therefore we add an energy window cut. In fig. 14, we see the sensitivity on BR for exotic Z decay for topology 6A can reach ∼ 10 −7 for Giga Z and 3 × 10 −9 for Tera Z. For m φ < 2 GeV, it is hard to separate the two photons from φ decay and the signal efficiency goes to zero. Instead of three photons in the final state, one could look for two photons because the photons from m φ can not be distinguished and therefore cover this mass range as in fig. 8.

IV.8. The sensitivity reach of the HL-LHC
The HL-LHC (3 ab −1 ) also produce a lot of Zs, which can be sensitive to some of the exotic Z decay modes. In this section we would like to study the sensitivities for the HL-LHC, and compare its reaches with the ones from Z-factories. A full fledged study with realistic detector simulations is beyond the scope of this paper. Instead, we perform simplified simulations aiming at gain an order of magnitude estimation. As we will see, the capabilities of the HL-LHC and Z-factories are very different. Our approach is sufficient to highlight the relative strengths of the two experiments. For each topology we only pick up one benchmark parameter (zero for DM mass and 40 GeV for other new physics mediate particles) to set the HL-LHC sensitivity, and did not scan the parameter spaces of models, because the cut efficiency is not strongly depending on the mass. We have chosen the benchmark mass parameters to give the most energetic Z-decay products. In addition, we do not consider fake photons from QCD, which will significantly reduce the HL-LHC sensitivity. In this sense, our projection for the HL-LHC should be considered as optimistic. In order to suppress the huge QCD background and avoid pre-scaling, we search Z production in association with a high p T jet or high p T photon. For all the visible particle, we require |η| < 2.5. Z → γ + / E: For exotic decay Z → γ + / E, we generate jZ event with the Magdraph 5 at 13 TeV LHC, and require the jet to be p j 1 T > 60 GeV to make Z have enough p T produce the energetic photon and large enough / E suppress the SM background. Specifically, we require / E T > 50 GeV and p γ T > 20 GeV together with p j 1 T > 60 GeV as the basic cuts. After the parton level event generation, it is passed to Pythia v6.4 [133] for showering and hadronization, and to Delphes v3.2 [134] for detector simulation. In the detector, missing energy could come from the jet reconstruction due to jet energy resolution and uncertainty. Therefore, we include the SM background jγ and irreducible SM background jγνν. We list the cross-sections after basic cuts for signal jZ → j + γ + / E and each SM backgrounds in table II in the column labeled with "σ basic ".  Table II. The exotic Z decay final states are listed for both SM backgrounds and signals. The "σ basic " column gives the cross-section after basic cuts, and the gives the cut efficiency for the further optimized cuts. The above cut efficiencies do not including the b-tagging efficiency. In the final sensitivity calculation, we use the b-tagging efficiency 0.7 and mis-tagged efficiency 0.015 [135] to re-weight the events according to the signal. For SM background jγγ, jγγνν, j + − and j + − νν, there is an additional invariant mass window cut for γγ or + − , which should multiply the efficiency given in the parentheses (×()). This additional efficiency is given as a range, because the mass window changes with the mediator mass in the signal topology. Such change is indicated by the light brown shaded region for HL-LHC in fig. 16.
To further optimize the signal, we make the differential distribution for kinetic variables p j 1 T , / E T and p γ T in fig. 15. We compare the distribution of SM background jγ and jγνν with signal jZ with exotic Z decay topology 1A, 1B and 1C. Based on fig. 15, we further impose the following cuts, We did not use additional cuts on p j 1 T because the distributions of SM background and signal are quite similar. After applying the above cuts, we list the corresponding cut efficiency in table II in the column labeled " ". For the HL-LHC (3 ab −1 ), we can reach the sensitivity for exotic Z decay BR of 5.6 × 10 −6 , 2.3 × 10 −5 and 5.76 × 10 −6 for signal topology 1A, 1B and 1C. The sensitivities for the HL-LHC for each topology are given in the summary plot fig. 16 Figure 15. The normalized event distributions for kinematic variable p j1 T , p γ T and / E T for signal jZ → j + γ + / E and the corresponding SM background. The distributions have been normalized to 1.
The SM background we consider are jγγ with / E from mis-reconstruction and irreducible jγγνν. The basic cuts are p j T > 60 GeV, / E T > 50 GeV, and two photons with p γ T > 20 GeV. The crosssections after cuts for signal and SM background are again listed in table II. We use the following cuts to further optimize our signal, The cut efficiencies for signal and SM background are listed in table II. For topology 2A and 2C, we can make an additional 5 GeV window cut on the invariant mass of diphoton to suppress the SM background, while the signal is remain unaffected. The corresponding efficiency has been listed in the parentheses in the column in table II. It is a range for SM background due to the change of mediator mass. For the HL-LHC (3 ab −1 ), the future sensitivity reach for exotic Z decay topologies 2A, 2C and 2D are (5 − 10) × 10 −7 , (1 − 2) × 10 −6 , and 1.4 × 10 −6 respectively, and have been plotted in fig. 16. The sensitivity range for the topology 2A and 2C has been indicated by light brown shaded region. Z → + − + / E: For decay topology Z → + − + / E, we consider SM background j + − and irreducible j − + νν with the same reason. The basic cuts are one jet with p j T > 60 GeV, missing energy / E T > 50 GeV, and two leptons with p T > 20 GeV. After checking the kinematic variable distribution, we propose further cuts to optimize our signal, For topology 3A and 3B, we have added the same additional 5 GeV window cut on the invariant mass of dilepton. The corresponding efficiency has been listed in the parentheses in the column in table II. For the HL-LHC (3 ab −1 ), the future sensitivity reach for exotic Z decay topologies 3A, 3B, 3D and 3F are (3 − 11) × 10 −6 , (3 ∼ 12) × 10 −6 , 2.0 × 10 −5 and 1.6 × 10 −5 respectively, and have been plotted in fig. 16. The sensitivity range for the topology 3A and 3B has been indicated by light brown shaded region in fig. 16. Z → jj + / E: For decay topology Z → jj + / E, we generate signal events γZ to suppress QCD background and consider SM background γj and irreducible γjjνν. The basic cuts are two jets with p j T > 30 GeV, missing energy / E T > 50 GeV, and one photon with p γ T > 60 GeV. After checking the kinematic variable distribution, we propose further cuts to optimize our signal, For the HL-LHC (3 ab −1 ), the future sensitivity reach for exotic Z decay topologies 4A, 4B, 4C are 0.0136, 3.45 × 10 −3 and 5.07 × 10 −3 , respectively, and have been plotted in fig. 16. Z → (JJ)(JJ): For decay topology Z → (JJ)(JJ) which is fully hadronic, we generate signal events γZ to suppress QCD background and consider SM background γJ matched with γJJ by Pythia and irreducible γJνν matched with γJJνν, where J can be light flavor jets j or b-tagged jet b. We require at least four jets with p J T > 60 GeV, and one photon with p γ T > 60 GeV. We propose further cuts to optimize our signal, and the cut efficiencies for signal and SM background are given in table II. Note we have generated the SM backgrounds with light flavor jet and b-jet separately. Both of them can contribute to background of the corresponding signal topologies 5A, 5B and 5C with b-tagging efficiency reweighting. For the HL-LHC (3 ab −1 ), the sensitivity reach for exotic Z decay topologies 5A, 5B, 5C are 0.0126, 0.0172 and 0.00915, respectively, and have been plotted in fig. 16. It is not surprising that sensitivity for fully hadronic decay of Z at the HL-LHC can not compete with future e + e − collider, because of the huge QCD background. Using jet substructure technique can probably achieve better sensitivities in the exotic hadronic Z decay topologies. CMS at 13 TeV LHC has searched for light vector resonance which decay into quark pair in association with a high p T jet to make the light vector gauge boson highly boosted [136], which decay products are merged into a single jet. The characteristic feature of the signal is a single massive jet with two-prong substructure produced in association with a jet from initial state radiation. The SM process jZ → j(jj) has been nicely reconstructed. In the exotic decay topology Z → (JJ)(JJ), one would look for four-prong substructure in the fat jet to suppress the SM jZ background. For final state including b-jets, b-tagging techniques in the jet substructure could further help in reducing the SM QCD background, which already help observing with a local significance of 5.1 standard deviations for the first time in the single jet topology in Z → bb process [137]. Z → γγγ: The last exotic Z decay search is Z → γγγ, which has been performed by ATLAS at 8 TeV LHC [122] with L = 20 fb −1 . The corresponding constraint is BR(Z → γγγ) < 2.2 × 10 −6 . It is hard for us to reliably study this topology due the difficulty in simulating the fake photons from QCD backgrounds. Instead, we do a simple rescaling according to the HL-LHC integrated luminosity 3 ab −1 , which gives the limit BR(Z → γγγ) < 1.8 × 10 −7 . Summary: In fig. 16, we see the sensitivity of exotic Z decay branching ratios at the HL-LHC generally can not compete with future Z-factory because of the large QCD background. The Z exotic decay products are typically rather soft. Requiring another hard radiation can help with triggering and making the Z-decay products more energetic. At the same time, it will reduce the signal rate significantly. For exotic Z decay with missing energy in final state, another important background can come from mis-measurement of the QCD jets. Since the missing energy from Z-decay tends to be small, this background can be significant. For photons in final states, there can be fake photons from QCD which we have not considered. For hadronic exotic Z decay, the situation at the HL-LHC is even worse.  Figure 16. The sensitivity reach for BR for various exotic Z decay topologies at the future Z-factory (Giga Z and Tera Z) and the HL-LHC at 13 TeV with L = 3 ab −1 . The BR sensitivity generally depends on model parameter, for example mediator mass and dark matter mass. The dark color region with solid line as boundary indicates the worst reach for the topology, while the lighter region with dashed line indicates the best reach. For HL-LHC, we add the light shaded region for the topology 2A, 2C , 3A and 3B to indicate the effect of an invariant mass window cut for diphoton and dilepton. For the topology 6A, the HL-LHC limit is obtained by rescaling the ATLAS study at 8 TeV LHC [122] with L = 20 fb −1 .

V. CONCLUSION AND DISCUSSION
We have presented a comprehensive study on exotic Z decay at future Z-factories, with emphasis on its prospects to exploring dark sector models. There are many dark sector models can give rise to exotic Z decay modes, many of which contain missing energy in the final states. A Z-factory provides a clean environment for decay modes which can be overwhelmed by large background at hadron colliders. Another advantage of searching for such exotic Z decay at future e + e − colliders is the ability of reconstructing the full missing 4-momentum, while we can only reconstruct missing transverse momentum at hadron colliders. We have demonstrated the capability of exotic Z decay at future Z-factory to provide the leading constraint in comparison with existing collider limits, future HL-LHC projections, and current DM searches.
We classify final states of the exotic decays with the number of resonances, and possible topologies it could have. We make projections on the sensitivity on the branching ratio of exotic Z decay at future Z-factory. For final states with missing energy, it can provide limits on BR down to 10 −6 − 10 −8.5 for Giga Z and 10 −7.5 − 10 −11 for Tera Z. The sensitivities on BR for different final states are roughly ordered from high to low as / E + − ∼ / Eγγ, / EJJ and / Eγ, due to the size of the SM backgrounds for each mode. In the same final states, it is quite clear the SM backgrounds for signal with more resonances can be better suppressed. In addition to the final states with missing energy, we also selectively studied the fully visible final states (JJ)(JJ) and (γγ)γ, where the first one contains two resonances and the second one contains one resonance. It is interesting to look for purely hadronic final states at future Z-factories, because it has much less QCD background in comparison with hadron collider. We found it can provide limits on BR down to 10 −5 − 10 −6.5 for Giga Z and 10 −6.5 − 10 −8 for Tera Z. The sensitivity to the final states with b jet is better than those with light flavor jets due to smaller SM backgrounds.
We have also made estimates of the reach of the HL-LHC on the exotic Z decay modes. The decay products tend to be soft and difficult for the LHC searches. There is also large QCD backgrounds. We considered the cases with additional energetic initial state radiation, which can help with suppressing these backgrounds. However, this also reduced the signal rate. Therefore, for the channels we considered here, it is very hard for the HL-LHC to compete with future Z-factory. The one exception is the (γγ)γ channel. The sensitivity on BR can reach 10 −7 for Giga Z and a few 10 −9 for Tera Z. The corresponding HL-LHC sensitivity is rescaled from an exiting study at 8 TeV by rescaling, and it can be comparable to that of the Z-factory.
We have studied four representative models in section III, namely Higgs portal with scalar DM, vector portal DM model, MIDM and RayDM, and axion-like particle model. In Higgs portal model with DM, the decay topology Z →sZ * → (χχ) + + − has been studied. Future Z-factories can provide the leading constraint on mixing angle sin α between SM Higgs and dark singlet scalar mediator. The constraint from Z →sγ via loop effect has also been considered, but is weaker due to loop suppression and larger SM background. In vector portal DM model, the decay topologies Z →Ã SS * → ( + − ) / E andZ →Ã φ → + − ( / E) are studied. The first one simply arises when DM is a scalar and charged under U (1) D , and the second one is a dark Higgs bremsstrahlung process. We found that the limits from the exotic Z decay provides a competitive and complementary constraints with DM direct detection, while the other collider limits are much weaker. In MIDM and RayDM model, the decay topologies Z → χ 2 χ 1 → (χ 1 γ)χ 1 from MIDM operator and Z → χ 1 χ 1 γ from RayDM operator has been considered. Both operators can be originated from heavy fermions and scalars in the loop, which couples to DM. The constraint on MIDM operator is much stronger than the constraint on RayDM. It is also much better than gamma-line search in indirect detection and future hadron collider projections. In axion-like particle model, the decay topologies Z → aγ → (γγ)γ and Z → aγ → ( / E)γ have been considered, where in the first one the axion-like particle decay promptly into two photons and in the second one it decays outside the detector. We find future Z-factory can provide the leading constraint on Λ aBB comparing with limits from LEP and LHC.
All in all, the exotic Z decay searches can provide unique tests on dark sector models at future Z-factory, especially when missing energy and/or hadronic objects appears in the final states. We explicitly analyze four representative dark sector models and find the exotic Z decay searches can provide the leading and complementary limits to the current and future collider searches and DM searches. It can also cover parameter spaces of DM models with the relic abundance requirement, which provides a complementary cross-check on DM problem.

VI.1. The annihilation cross-section for scalar DM with vector portal
We calculate the annihilation cross-sections of scalar DM into SM fermions. The scalar DM is charged under U (1) D as in eq. (17), and the kinetic mixing induced interactions with SM sector are given in eq. (16), which includes both s-channelÃ andZ mediation. The annihilation crosssections for one generation are given, where we see the cross-sections are p-wave suppressed.