Searching for lepton portal dark matter with colliders and gravitational waves

We study the lepton portal dark matter (DM) model in which the relic abundance is determined by the portal coupling among the Majorana fermion DM candidate $\chi$, the singlet charged scalar mediator $S^\pm$ and the Standard Model (SM) right-handed lepton. The direct and indirect searches are not sensitive to this model. This article studies the lepton portal coupling as well as the scalar portal coupling (between $S^\pm$ and SM Higgs boson), as the latter is generally allowed in the Lagrangian. The inclusion of scalar portal coupling not only significantly enhances the LHC reach via the $gg\to h^*\to S^+S^-$ process, but also provides a few novel signal channels, such as the exotic decays and coupling deviations of the Higgs boson, offering new opportunities to probe the model. In addition, we also study the Drell-Yan production of $S^+S^-$ at future lepton colliders, and find out that the scenario where one $S^\pm$ is off-shell can be used to measure the lepton portal coupling directly. In particular, we are interested in the possibility that the scalar potential triggers a first-order phase transition and hence provides the stochastic gravitational wave (GW) signals. In this case, the terrestrial collider experiments and space-based GW detectors serve as complementary approaches to probe the model.


Introduction
The Standard Model (SM) of particle physics has been a great triumph in explaining and predicting the astrophysical and terrestrial experimental phenomena, however there are still many unsolved problems remaining, such as the dark matter (DM). Many astrophysical evidences support the existence of DM, and the fitting result of Cosmic Microwave Background (CMB) to the ΛCDM model yields a DM relic abundance of Ω DM h 2 ≈ 0.12 [1], which accounts for ∼ 27% of the total universe energy. However, we still know very little about the particle origin of DM, except that none of the SM particles can be the DM candidate [2]. Therefore, the existence of DM is a clear evidence for physics beyond the SM (BSM).
Over the past several decades, the most popular particle explanation for DM has been the freeze-out mechanism of the weakly interacting massive particles (WIMPs) [3], as it naturally yields the observed DM relic density when the coupling of DM to the SM particles is of the order of the electroweak (EW) gauge couplings, and the DM mass is O(100 GeV).
The lepton portal DM model is proposed in Ref. [23], which assumes a portal coupling among the SM lepton and two dark sector particles S (scalar) and χ (fermion), where "dark sector" means an odd Z 2 symmetry is assigned to S and χ. Depending on the mass hierarchy, the DM candidate could be S or χ, while the fermion DM case can be further classified into Dirac or Majorana DM scenarios. This particular model has been further studied in Refs. . In this article, we consider the model with right-handed lepton portal, taking χ as Majorana DM candidate and S as the charged scalar mediator, which is similar to the setup in Ref. [23].
This model has very small indirect search cross sections due to helicity suppression. Because the DM candidate only couples to charged leptons, its nuclear recoil cross section comes from loop diagrams. Thus, the direct search signal is also suppressed. Therefore, collider experiments are crucial in probing this model. The typical collider signal is the Drell-Yan pair production of the S ± mediator, and its subsequent decay to χ and a charged lepton. In this article, we study this channel at the LHC and future e + e − colliders, and in the latter case we include the off-shell S pair production S ± S ∓( * ) , which provides a direct probe for the lepton portal coupling.
Different from previous studies, in addition to the lepton portal coupling, we also consider the Higgs portal coupling |S| 2 |H| 2 , which is in general allowed in the Lagrangian. The inclusion of this coupling leads to several novel signals, such as the gluon-gluon fusion production of S + S − at the LHC, the Higgs exotic decay (e.g. h → χχ, h → + − χχ), the Higgs coupling (e.g. hZZ, hγγ, h + − ) deviations and the lepton (g − 2) corrections to the SM prediction. In particular, the scalar portal coupling might be able to trigger a firstorder phase transition (FOPT) in the early universe, opening the window for detecting the model via the stochastic gravitational wave (GW) signals. We finally consider the interplay between the GW searches and collider searches and show their complementarity. This paper is organized as follows. We describe the model and derive the parameter space for the WIMP DM candidate in section 2. In section 3, various terrestrial searches are investigated, including the S + S − production, the exotic decay and coupling deviations of the Higgs boson at the LHC and future e + e − colliders, together with lepton (g − 2) searches to constrain the two portal couplings in this model. The scenario that this model provides the right DM relic abundance and at the same time triggers a FOPT is considered in section 4, in which we show that detectable GW signals suggest a large Higgs portal coupling. In section 5, the interplay between GW detectors and the collider experiments are discussed. Finally, we summarize and conclude in section 6.  Figure 1. The lepton portal coupling y th as a function of m S and m χ , which satisfies the DM relic abundance requirement.

The model
The model contains two new fields: the Majorana DM candidate χ, which is a gauge singlet; and the complex scalar mediator S, which is an SU (2) L singlet with hypercharge −1. The relevant Lagrangian reads where H is the SM Higgs doublet, and = e, µ, τ is the SM charged lepton (mass eigenstate). We require S couple to one generation of lepton at a time to avoid lepton flavor violation. Such a flavor alignment typically needs some specific underlying mechanism to realize, which however beyond the scope of this work. For a more generic model that simultaneously involves all three flavors, our results still apply, as long as suitable rescaling is performed. The model contains a Z 2 symmetry for χ and S, that both of them carry odd charges. Assuming m χ < m S , the χ is stable and thus can be the DM candidate.
The χ pair can annihilate into the lepton pair via the exchange of a t-channel S. Due to the Majorana nature of χ, the s-wave component of χχ → + − is suppressed by the lepton mass. Therefore, the annihilation cross section is p-wave dominated [65,66], where x ≡ m 2 χ /m 2 S , and we have applied the limit m → 0 in the second equality. We perform the thermal average of the annihilation cross section according to Ref. [67]. Given a set of (m S , m χ ), one can always tune y to have the right annihilation cross section at freeze-out to achieve the observed DM relic abundance, and the corresponding y is denoted as y th , which is plotted in Fig. 1. One can see that for EW scale m χ and m S , a Yukawa coupling y ∼ O(1) can provide the correct DM relic density. For m χ < m S , smaller m χ leads to larger y th , because the annihilation cross section scales as m 2 χ /m 4 S . As the DM annihilation signal χχ → + − is helicity or p-wave suppressed, it is hard to be probed by satellite experiments like Fermi-LAT [68,69], AMS-02 [70,71], or the CMB measurements from Planck [1]. For direct detection, the scattering between χ and nucleons arises only at one-loop level, which can be described by an effective operator. Since χ is a Majorana fermion, its dimension-five magnetic dipole operator vanishes. It leaves the dimension-six operator as the leading contribution, which can be matched to the electromagnetic anapole moment of DM. This receives additional suppression from DM velocity square, so that it is difficult to detect from the direct detection experiments [23]. 1 As a result, we conclude that the lepton portal DM with Majorana DM has negligible signal in indirect and direct searches, and is only subject to the constraints from the thermal relic abundance and collider searches.

The particle experiment searches
In this section, we discuss probing the lepton portal coupling y and Higgs portal coupling λ HS in particle experiments. First, we consider the S + S − → + χ − χ channel at the LHC and future lepton colliders. In the latter case, the off-shell production of S ± offers the opportunity to probe y directly. Next, we study the exotic decays of the Higgs and Z bosons, including the three/four-body decays h/Z → S ±( * ) S ∓( * ) and the invisible decay h → χχ at the Higgs factory CEPC and FCC-ee, which probe the combination of couplings y and λ HS . Then we turn to the correction to the Higgs couplings h + − , hγγ and hZZ. Finally, the lepton (g − 2) is discussed.
3.1 Pair production of S ± 3.1.1 pp → S + S − at the LHC In the model, the lepton portal scalar S carries one unit of hypercharge, therefore it can be produced in pair via the EW Drell-Yan process qq → Z * /γ * → S + S − mediated by off-shell γ and Z bosons. However, in our model due to the large scalar sector coupling λ HS , one can also have the (off-shell) Higgs mediated process pp → h * → S + S − which can significantly modify the total cross section of pp → S + S − . For example, with λ HS = 1 and m S = 200 GeV, the cross section contributed by the Higgs mediation can be 30% of the total cross section. The production rates of the Drell-Yan and the gluon-gluon fusion processes are shown in Fig. 2.
The produced S ± exclusively decays to ± χ, leading to a di-lepton plus missing transverse energy final state ( + − + / E T ) at the LHC. The LHC experiments have already set constraints on such a final state in searches for the sleptons from SUSY models [73][74][75]. The LHC Run-I and Run-II data from ATLAS [73,74] have covered mass of S ± up to 450 GeV for the exclusive decay channel e ± χ or µ ± χ. The compressed parameter region when 1 The low energy electron recoil cross section is (y 4 /π)m 2 e /(m 2 φ − m 2 χ ) 2 , which is typically 10 −44 cm 2 for m φ , mχ ∼ O(100) GeV, well below the constraint from LUX-ZEPLIN experiment [72]. m χ is close to m S , has also been studied by the ATLAS collaboration [75]. Earlier studies from LEP have fully excludes such charged scalar S with mass m S < 100 GeV [76]. We show the above existing constraints in Fig. 3 as colored regions.
To make future projections for our model, we write down the UFO model file [77] with the FeynRules package [78] and use MadGraph5 aMC@NLO [79] to generate parton level events, and then use Pythia8 [80] and Delphes [81] to implement parton shower and fast detector simulation, respectively. Both the EW Drell-Yan and the Higgs mediated S + S − production processes are included in our simulation. The signal events are selected using the cuts from two signal regions (SRs) of the ATLAS study [74] as follows, Here the m T 2 observable is defined event-by-event as the minimum of the function [82] max m T ( p  According to the number (1 or 2) of light-flavor jets in the final state and the different m T 2 cuts (100 or 120 GeV), we classify the events into 4 SRs, and adopt the simulated background event numbers from Ref. [74]. We use S/ √ B = 2 as the criterion of the LHC reach for a given integrated luminosity, where S and B denote the signal and background event numbers, respectively. The projections at the 14 TeV LHC with 300 fb −1 are plotted in Fig. 3 as black lines, where we show both the pure Drell-Yan contribution (means λ HS = 0, dashed) and the inclusion of gluon-gluon fusion results for λ HS = 3 (dotted). For large λ HS , the reach can be visibly enhanced. Note that the enhancement is not significant in the low m S region, although in that region the gluon-gluon fusion process has a larger fraction, as shown in Fig. 2. That is because the S ± from gluon-gluon fusion typically have a softer p T , and hence they are cut away by the hard m T 2 cut in our simulation. Loosing m T 2 might help to probe the low m S region, and we leave the detailed study for a future work.
3.1.2 e + e − → S ± S ∓( * ) at future e + e − colliders As shown in Fig. 3, there is a gap between the LHC and LEP constraints for 100 GeV < m S 150 GeV and 30 GeV m χ 100 GeV. The future e + e − colliders with a collision energy of ∼ 250 GeV can fill this gap. Moreover, an e + e − machine is able to probe the lepton portal coupling y directly, provided one S ± is off-shell. For the on-shell production at LHC, since S ± → ± χ decay branching ratio is 100%, the rate does not depend on y . Therefore, the exclusion of slepton-like particle S at LHC is shown only in the m S -m χ plot, but can not constrain y . However, for the 2 → 3 process e + e − → S ± ∓ χ mediated by an off-shell S ∓ , the rate does depend on y 2 , opening the window to directly probe the DM  The analysis is carried out on the typical Higgs factory such as CEPC, ILC or FCC-ee, with √ s = 250 GeV. We include the on-shell 2 → 2 S + S − pair production as well as the off-shell 2 → 3 process. For = e, the 2 → 2 process includes both Drell-Yan contribution and the t-channel χ mediated contribution. Therefore, in this special case, the 2 → 2 cross section already has the contribution of y e . When = µ, τ , the cross section does not depend on y , since only Drell-Yan process contributes. Regarding 2 → 3 process, there are two advantages to include it. One is probing the y coupling, and the other is that we can probe √ s/2 < m S < √ s region. Note that given m S and m χ , y is determined by the relic abundance of the DM.
We generate signal events using the packages mentioned in section 3.1.1. Following Ref. [83], we add a few cuts to select the signal events, i.e.
1. Exactly two opposite charged leptons with p T > 5 GeV and |η| < 3; 2. Veto any events with a jet within p T > 5 GeV and |η| < 3; The background event numbers are adopted from the simulated results from Ref. [83]. In Fig. 4, we can see that this search is complementary with the LEP and LHC results. It can cover the region 100 GeV < m S 150 GeV region, which is not touched previously when m χ mass is moderately large, e.g. 30 ∼ 100 GeV. For large m χ , the visible energy shared in the leptons decreases, which makes it hard to compete the LHC background from W + W − pair. The cleaner environment at the lepton collider makes it sensitive to softer leptons comparing with LHC. Including 2 → 3 process, we do see the sensitivity at future e + e − collider extends to m S ∼ 170 GeV and 150 GeV for y e and y µ respectively. For small m χ there is higher reach for m S , because in this region a large y is needed to get the correct DM abundance (see Eq. (2.4)), which enhances the signal significance at the collider (through the off-shell S ± ).

Exotic decays from the Higgs and Z bosons
A charged S ± with m S < 100 GeV is already excluded by the LEP experiment, forbidding the h → S + S − one-shell decay. However, for m χ < m h /2 ≈ 62.5 GeV, the λ HS portal coupling can induce the exotic three-or four-body decays h → S ± ∓ χ or h → ± χ ∓ χ mediated by one or two off-shell S ± , depending on whether m S < m h or not. The decay width is proportional to y 2 λ 2 HS or y 4 λ 2 HS , providing a new way to probe λ HS and y . If we fix y = y th by the relic abundance requirement, it gives a limit to λ HS for a given set of (m S , m χ ).
We explore this exotic Higgs decay at the Higgs factory FCC-ee and CEPC at √ s = 250 GeV via the e + e − → Zh production channel, whose cross section is 0.24 pb. We consider the following cascade decays Z → + − and h → S ±( * ) S ∓( * ) → + χ − χ, 2 where χ plays the role of missing energy. The main SM backgrounds include ZW + W − and Zτ + τ − , with all the particles decaying to leptonic final states, and the Zh, with h decaying to W ± W ∓ * , ZZ * and τ + τ − . The total cross section for the backgrounds in the ( + − ) + − / E T final state is as small as 0.02 fb, where the leptons in the parentheses come from Z decay. The signal and background events are simulated by the packages mentioned in section 3.1.1. We apply the following detailed requirements to the events, 1. At least four charged leptons with p T > 10 GeV and |η| < 2.47; 2. A pair of leptons with same-flavor and opposite-sign, and satisfies |m + − − m Z | < 5 GeV; 3. The missing energy / E T > 20 GeV; 4. The reconstructed Higgs resonance in the mass window [84] 120 which is equivalent to cut on the total energy of Z.
After the cuts, for a given integrated luminosity we are able to set bounds for Br(h → S ±( * ) S ∓( * ) → + χ − χ), which in turn is translated into the upper limits for λ HS once y is fixed by the relic abundance requirement, as shown in left panel of The lepton portal coupling y is set to be y th , which is the value to provide the correct DM relic abundance. Right: The 95% C.L. constraint contours (magenta) for y for the exotic decay Z → + χ − χ. The gray shaded region shows the parameter space can be excluded if y = y th .
of the curves around m S = 125 GeV is originated from the phase space change from threebody to four-body decay. In conclusion, the future e + e − collider can significantly constrain the couplings for light m χ 30 GeV and m S of a few hundreds GeV. In addition to the exotic SM Higgs decay, another good target to probe new physics is the Z exotic decay [85]. To explore the SM parameters with better precision, the future e + e − colliders have the proposals to run at Z-pole [86][87][88], which can provide Giga (10 9 ) or Tera (10 12 ) Z boson. Interestingly, in this model, there exists the exotic decay channel Z → + χ − χ, which provides a di-lepton plus missing energy final state. There are two types of diagrams responsible for this channel, the first one involves two off-shell S ± through ZS + S − coupling and the second one involves one off-shell S ± through Z + − coupling with S attached to one of the charged lepton. Since we are considering S ± with mass larger than 100 GeV, it is heavier than all other particles. Therefore, the decay width is dominated by the latter diagram, which is proportional to y 4 m −4 S . We explore this exotic Z decay similar to the Higgs case. The dominant SM background isνν + − from off-shell gauge boson pair production. 3 The cut conditions are: 1. At least two same-flavor opposite-sign leptons with p T > 10 GeV and η < 2.5; 2. The missing energy should satisfy / E T > 10 GeV.
The 95% C.L. constraint on the exotic decay branching ratio is about Br(Z → e − e + χχ) 10 −9 for Tera Z option, while exact limit again depends on m χ and m S . In right panel of shows that for Tera Z option, the thermal DM with mass m χ 13 GeV can be excluded by this exotic Z search, because y th exceeds the limit from the exotic Z decay Z → e − e + χχ.
We plot such region in gray and label as "thermal DM excluded". This constraint provides a complimentary limit for large m S comparing with LHC limits, because it does not require on-shell S production. Moreover, since the decay width and the DM annihilation cross section are proportional to y 4 m −4 S , this exclusion line can extend horizontally to very high m S , thus is a powerful constraint for this DM model.

Invisible decay: h → χχ
The Higgs invisible decay h → χχ is induced at one-loop level by the two Feynman diagrams listed in Fig. 6. 4 The first diagram is negligible due to the small lepton mass. We calculate the second diagram contribution to the exotic Higgs decay using Package-X [92]. The induced coupling is where DiscB is the finite part of the Passarino-Veltman function B 0 defined in Package-X, C 0 is the Passarino-Veltman function following the definition in Package-X and the lepton mass is taken to be zero. In the second equality, we have expanded the result with large m S . We have checked our result with Ref. [90] and Ref. [35], and found agreement between each other. In numeric calculation, we use the full expression from Eq. (3.5).
The h → χχ partial width is given as (3.6) The current best limit for invisible Higgs decay is Br(h → inv) < 13% at ATLAS Run-II with integrated luminosity 139 fb −1 [93]. For future HL-LHC, the sensitivity for invisible Higgs BR is 3.5% from [94]. For future e + e − collider such as CEPC, the sensitivity can be increased to about 0.3% [86]. The above projections from future colliders can set limits on the coupling combination y 2 λ HS as a function of masses m S and m χ . In the left panel of Fig. 7, we show the contours of y 2 λ HS for LHC (brown), HL-LHC (blue) and CEPC (red) sensitivities. The dashed and solid contours corresponds to y 2 λ HS = 1, 10 respectively. We clearly see that the future e + e − collider has a better sensitivity over the hadron colliders. In the right panel of Fig. 7, y is set to be y th to satisfy DM relic abundance requirement. Once fixing m S and m χ , the future sensitivity on λ HS can be calculated using CEPC sensitivity Br(h → inv) = 0.3% and the contours are shown. One interesting feature is that for m χ smaller than 6 GeV, the sensitivity on λ HS goes down significantly, because the leading order in the width for small m χ expansion is linear in m χ . When m χ decreases further, y th increases to compensate the annihilation cross section, which makes the sensitivity on λ HS becoming stronger again. Therefore, we see λ HS reaches its best sensitivity for small m S and moderate m χ .

One-loop contributions to Higgs couplings
The h + − vertex is modified by the one-loop diagram in Fig. 8, 5 and the induced interaction is calculated using the Package-X where we have taken the leading power in small lepton mass m . It is obvious that the coupling is lepton mass suppressed, because the DM only couples to right-handed lepton 5 Similar corrections have been studied in SUSY model already [90,95]. Figure 8. The one-loop induced h + − coupling modification. branching ratio and coupling strength κ τ from existing LHC constraints [97,98] with fixed DM mass m χ = 10 GeV. Middle and right panels: the constraint on λ HS from the projected precision δκ µ < 8.7% and δκ τ < 1.5% from CEPC [99], with y = y th by DM relic abundance. and one has to flip the helicity of lepton, which induce this suppression. Of course, the SM Higgs couplings to leptons are suppressed by the mass as well. Thus, the one-loop correction to the SM coupling in fraction is proportional to y 2 λ HS /(8π 2 ). In addition to vertex correction, the Higgs field renormalization and the lepton field renormalization can also contribute to the coupling. The experiments characterize the sensitivity of Higgs coupling measurements in the κ-framework [96], which is calculated using the cross section σ(Zh) and H decay branching ratios and hence the universal Higgs field renormalization effect cancels out in the decay branching ratios. Therefore, Higgs field renormalization does not affect the branching ratio to leptons. For the lepton field renormalization, it only contributes to right-handed lepton. We calculated that it is about y 2 /(128π 2 ) for large m S expansion, which is much smaller than Eq. (3.7). Together with the fact that λ HS 1, we can neglect the right-handed lepton field renormalization effect and focus on the vertex correction.
The results are shown in Fig. 9. In the left panel, we show the limits for the combination y 2 λ HS from Br(h → µ + µ − , τ + τ − ) at LHC and CEPC. The shaded regions are already excluded by current LHC limits, while the future sensitivities from CEPC are plotted as dotted lines. The DM mass is fixed as m χ = 10 GeV and the limits get weakened linearly with increasing m S . In the middel and right panels, we fix y = y th for the DM relic abundance requirement, and show the sensitivity contours for λ HS as a function of m S and m χ . When taking m S → ∞, the coupling in Eq. (3.7) is proportional to y 2 m −2 S . Therefore, it has the similar dependence as the DM relic abundance requirement which is y 4 m −4 S . This feature is clearly shown that when m S increases, the contours for λ HS are flat. It means that the sensitivity does not suffer for large m S , because the large mass m S is compensated by large y th . Since for µ and τ , this change of coupling ratio δκ does not depend on lepton mass, the constraint on λ HS is only linear depends on the CEPC precision, which has better τ sensitivity. Therefore, the sensitivity for τ lepton portal is better by a factor of 5.8 compared to µ lepton portal.
Besides the one-loop modification to h + − coupling, the charged scalar loop can also modify the hγγ and the hZZ couplings [100]. The former one modifies Br(h → γγ) via the charged scalar triangle loop, and current fitting result for the signal strength is 1.12 ± 0.09 at the 13 TeV LHC with 137 fb −1 [101], while the projection accuracy is 4% at the HL-LHC [96], better than the projected sensitivity 6.9% at the 5.6 ab −1 CEPC [86]. The latter one is mostly contributed by Higgs field renormalization (via the S ± loop) if λ HS > e (the elementary charge). This modification doesn't rely only on y , therefore it provides a direct constraint on the scalar interaction coupling λ HS independent of the Yukawa coupling. Current constraint for δκ Z is 8% for the 13 TeV 139 fb −1 LHC [102], and the projected result is 1.7% at the HL-LHC [96]. Therefore, the precision measurement on Zh production cross section can provide a limits on λ HS . Following the procedure of Ref. [100], we calculate the sensitivity for λ HS from the future limits on Br(h → γγ) and the coupling strength to Z gauge boson κ Z . The results are presented in Fig. 10.

The lepton g − 2
Another consequence of the lepton portal coupling y is the magnetic dipole moment for leptons at one-loop level. The g − 2 contribution for lepton is given as [103,104] where x ≡ m 2 χ /m 2 S and we keep the leading order result in the limit of small m . The last term containing x in Eq. (3.8) is a monotonically decreasing function of x and it equals to 1/6 and 1/12 in the limit of x → 0 and x → 1. Moreover, ∆a is always negative. The electron magnetic dipole moment has been directly measured in Refs. [105,106], and one can compare it with the SM prediction [107,108]  where ∆a e comes from the α measurement in recent cesium recoil experiments [109], while ∆ã e comes from a new independent measurement using rubidium atom [110]. For electron g − 2, the uncertainties from QED, EW and hadronic contributions are much smaller than α, while for muon g − 2, the uncertainty in α is less significant. ∆a e has a mild anomaly at confidence level of 2.4 σ, which might need a negative contribution from new physics. However, ∆ã e turns to be positive and the significance goes down to 1.6 σ. The measurement [110] placed 95% C.L. bounds to be ∆ã e ∈ [−34 × 10 −14 , 98 × 10 −14 ] and we will adopt this value in the later calculation. In our model, the contribution is always negative for the electron g − 2. In Fig. 11, We consider the y e range which can fit to ∆a e (blue shaded region) and ∆ã e (red shaded region) at 95% C.L., with a fixed m χ which is much smaller than m S . We fix m χ = 1 GeV for (g − 2) e fit, which satisfies the condition m χ m S easily. For g − 2 contribution in this limit, one can have the following expansion with small x as Therefore, for m χ considered in Fig. 11, we see that ∆a is determined by m S and y most of the time. As a result, the blue and red shaded regions will not shift when changing to other m χ . On the other hand, we plot y th e with a given DM mass, which can provide the correct DM relic abundance in Fig. 11. They are plotted as dot-dashed gray lines with a fixed m χ . We can see that for ∆a e , DM mass between 0.2 ∼ 2 GeV is preferred and can satisfy DM relic abundance at the same time, while for ∆ã e , DM mass larger than 1 GeV is required. The two measurements have a mutual region when DM mass is between 1 and 2 GeV.
The muon g − 2 also has a long-standing discrepancy [111,112]. Combination of the newest Fermilab and previous BNL measurements yields [113] ∆a µ = (2.51 ± 0.59) × 10 −9 , (3.11) corresponding to a significance of 4.2σ, suggesting a positive contribution from the new physics. 6 Unfortunately, in this lepton portal model the contribution has a negative sign that it can not explain the anomaly. To incorporate the ∆a µ result, one has to add new ingredients beyond current model.

Probing the model with gravitational waves
In this section we investigate the possibility of probing the scalar sector via the GWs from a FOPT in the early universe. 7 The first subsection is devoted to the discussion of FOPT, while the second subsection studies the GW detection limits.

First-order phase transition
The scalar potential in Eq. (2.1) is In terms of the real components, we get There are debates on this "excess". A lattice group shows that there is no significant tension between the SM prediction and the recent FNAL experimental determination [114]. 7 The FOPT of such a singlet charged scalar was also studied in Ref. [115]. Also see Refs. [116][117][118] for the FOPT triggered by a complex gauge singlet.
Here the quartic coefficients should satisfy (4.3) to ensure the potential is bounded below. Since S has electric charge number −1, it cannot develop a vacuum expectation value (VEV) at zero temperature. Hence the vacuum configuration at zero temperature is along the Higgs direction, i.e. ( h , φ ) = (v, 0). This means that the Higgs-relevant coefficients have been fixed by the collider measurements, where m h = 125 GeV and v = 246 GeV, leaving only three free parameters in the scalar potential. If µ 2 S is positive, then the vacuum configuration (v, 0) is trivially achieved, because along the φ direction there is no local minimum. However, if µ 2 S < 0, then the φ direction also has a local minimum w = −µ 2 S /λ S . In this case, to make sure (v, 0) is a minimum but not a saddle point, the coefficients must satisfy according to the Hessian matrix [119]. In case of λ S µ 2 H > λ HS µ 2 S , (0, w) is also a local minimum, and a further condition is required to ensure (v, 0) is the global minimum. The mass of S is given by m 2 S = µ 2 S + λ HS v 2 . Although the collider experiments have set a bound for m S (which also depends on m χ ), they cannot probe µ 2 S directly and neither the sign of µ 2 S . If µ 2 S < 0, the potential might trigger a FOPT in the early universe, and this opens a new detection scenario for this kind of DM models: the FOPT GWs.
At finite temperature, the scalar potential is modified by the thermal correction. Taking the leading gauge invariant T 2 terms [120,121], the thermal potential is where The necessary condition for a FOPT is the existence of a critical temperature T c at which the system has two energetically degenerate vacua (v c , 0) and (0, w c ). For Eq. (4.7), this requires [119] 9) and the critical temperature and VEVs are respectively Below T c , the Higgs-direction minimum becomes the lower one and the universe starts to decay to it from the φ-direction minimum. The decay rate per unit volume reads Γ(T ) ∼ T 4 e −S 3 /T , where S 3 is the classical action for the O(3) symmetric bounce solution [122]. The nucleation temperature T n is defined by the equality of the nucleation rate per Hubble volume and the universe expansion rate, i.e. Γ(T n ) = H 4 (T n ). For a radiation-dominated universe and a FOPT happening around the EW scale, T n can be solved by [123] S 3 /T n ∼ 140. (4.11) Eq. (4.11) is treated as the sufficient condition for a FOPT in this article. We use the Python package cosmoTransition [124] to calculate the bounce solution and hence T n for the potential in Eq. (4.7). As discussed before, after taking into account the Higgs mass and VEV, in the scalar potential there are only three free parameters, which we choose to be λ S , m S and λ HS . For a fixed λ S , the parameter space allowed by the FOPT can be projected to the m S -λ HS plane. We plot the parameter space for different λ S as shaded areas in Fig. 12. The shapes of those areas can be well understood by the analytical relation Eq. (4.9). A sizable λ HS is preferred by FOPT, because the phase transition is triggered by the potential barrier induced by the λ HS h 2 S 2 /2 term.

Gravitational waves
Stochastic GWs are produced during a FOPT via bubble collision [125,126], sound waves in the plasma [127] and the magneto-hydrodynamics turbulence [128,129]. The expanding bubble wall only accelerates for a short time before it reaches its final velocity v b < 1. Therefore, only a tiny fraction of FOPT energy deposits in the bubble shell, and the bubble collision contribution to the GWs is negligible. Instead, most of the phase transition energy is pumped into the surrounding fluid shells, making sound waves the dominant contribution to FOPT GWs [130]. However, the sound wave period only lasts for a finite time, after which the energy in the bulk fluid will cause the turbulence, which is another source for GWs [131]. Consequently, the GW spectrum today can be expressed as Ω GW (f ) = Ω sw (f ) + Ω turb (f ), (4.12) where f is the frequency, the subscripts "sw" and "turb" denote sound waves and turbulence respectively, and Ω GW is the ratio of GW energy density to the critical energy of the current universe, i.e.
where ρ c = 3H 2 0 /(8πG), with H 0 being the Hubble constant today. The sound wave and turbulence spectra can be expressed as functions of two FOPT parameters [132][133][134], where ∆V T denotes the (negative) effective potential difference between the true and false vacua, and g * ∼ 100 is the number of relativistic degrees of freedom during FOPT. 8 Namely, α is the ratio of FOPT latent heat to the radiation energy, while β/H is the inverse ratio of FOPT duration to the universe expansion time scale. Numerically, the sound waves spectrum is [127] Ω (4.16) To take into account the finite duration of the sound wave period, we make the replacement Ω sw (f ) → Ω sw (f )H(T n )τ sw in Eq. (4.15) [131,137], where A more accurate treatment for the sound wave cutoff factor can be found in Ref. [138].
The turbulence spectrum is [128,129] where 8 It is suggested that the α and β/H parameters should be calculated at the percolation temperature Tp [130,131,[135][136][137]. However, we have checked that for our FOPT scenario α 1, therefore the supercooling effect is not prominent and Tn ≈ Tp is a good approximation. (4.20) The factor κ v and κ turb in Eq. (4.15) and Eq. (4.18) are the fraction of FOPT energy that is transformed to bulk motion/turbulence, respectively. We have adopted κ turb = 0.05κ v , and κ v is extracted from the numerical function of Ref. [139]. v b = 0.6 is used as a benchmark in our study, although it might be calculated using the hydrodynamics. 9 The GWs from a FOPT around EW scale might be probed by the next generation space-based laser interferometers such as LISA [142], BBO [143], TianQin [144,145], Taiji [146,147] and DECIGO [148,149]. The detectability of GW signals can be characterized by the signal-to-noise ratio (SNR). Taking the LISA detector as an example, the SNR is defined as [133] where Ω LISA is the sensitive curve of the LISA detector [133] (in which we take the C1 configuration as a benchmark), and T = 9.46 × 10 7 s is the data-taking duration [134].
We adopt SNR = 10 as the detectable threshold, and found that for a fixed λ S , only a narrow band of parameter space in Fig. 12 can be probed. The parameter space that can be probed by LISA is plotted in the top panel of Fig. 13 where (0, w n ) and (v n , 0) are respectively the old and new vacua at the nucleation temperature T n .

The interplay between phase transition and particle searches
In this section, we revisit the parameter space for λ HS from GW study and the corresponding constraints from the collider studies. We focus on how the two different types of experiments can be complimentary with each other. For GW experiment, we focus on the LISA detectable parameter space in the left panel of Fig. 13. Clearly, the scalar S self-interaction coupling λ S can affect the LISA detectable parameter space. Therefore, we vary λ S between 0 and 4π to obtain the entire GW detectable region for a FOPT.

Interplay with pp → S + S − at the LHC
The S ± particles can be pair produced at the LHC via Drell-Yan and gluon-gluon fusion processes, and the current constrains and future projections at the 14 TeV LHC has been studied in Section 3.1.1. As shown in Fig. 3, a non-zero λ HS can visibly enhance the probe limit in parameter space due to the gg → h * → S + S − contribution to the signal events. On the other hand, a sizable λ HS may trigger a FOPT and hence give detectable GW signals as well.
The interplay between the LHC and the future LISA experiments is plotted in Fig. 14. The light blue (orange) shaded region corresponds to λ HS = 2 (3), with the vertical boxed boundary regions being the LISA-detectable parameter space, while the irregular boundary regions being the enhanced LHC projections when including the gg → h * → S + S − contribution. This is because the GW signals is independent of the DM mass  Figure 15. The interplay between gravitational wave detection and future e + e − collider searches. The gray shaded region is the LISA detectable parameter space, varying λ S from 0 to 4π. From left to right, we show the sensitivities for λ HS from future CEPC precision measurements, based on invisible Higgs decay branching ratio Br(h → inv) = 0.3%, Higgs leptonic coupling precision reaches δκ µ < 8.7% and δκ τ < 1.5%. m χ . For a given λ HS , there is a set of upper and lower bounds for m S in the LISAdetectable region. The enhanced parts of the LHC probe region due to λ HS are also shown in the figure. We see that the LHC and GW experiments mainly serve as complementary approaches to probe the DM parameter space; while they also have some intersections, which can be used to identify the origin of the excess (if found).

Higgs precision measurement at the future e + e − colliders
The collider constraints on λ HS have to be related with SM Higgs. The constraint from exotic Higgs decay is less sensitive compared to the Higgs 1-loop coupling as shown in the previous section. The 1-loop induced Higgs couplings include the coupling to χχ and + − . The former can be revealed by the Higgs invisible decay branching ratio, for example we consider the future sensitivity from CEPC Br(h → inv) = 0.3%. The latter can be revealed by the Higgs precision measurements at CEPC with relative precision of couplings δκ µ < 8.7% and δκ τ < 1.5%. In Fig. 15, we take the DM mass m χ = 1, 5, 10, 30, 50, 60 GeV respectively to show its effect on the sensitivities for λ HS . For a fixed DM mass, the corresponding colored line shows the maximum allowed λ HS from the future e + e − collider searches. In general, the exclusion power is better for light DM mass m χ .
In the left panel of Fig. 15, one can see that for m χ < 40 GeV the constraints on λ HS are quite similar. The reason is that the 1-loop induced coupling is proportional to y 2 λ HS m χ /m 2 S for large m S . At the same time, the annihilation cross section is proportional to y 4 m 2 χ /m 4 S which requires this combination to be a constant to satisfy the relic abundance. Therefore, the limits on λ HS from Higgs invisible decay is a constant. The colored lines in the left panel do show this feature, except when m S is too close to the Higgs mass and the expansion on large m S is not valid anymore, the sensitivity on λ HS changes slightly. For larger m χ , the sensitivity of λ HS is downgraded because the phase space suppression in the h → χχ decay. It is worth to mention that the sensitivity from Higgs invisible decay works equally good for all three flavors of lepton portals. This search can test most of the LISA detectable parameter regions for m χ < m h /2.
In the right panel of Fig. 15, the limits from δκ τ < 1.5% are plotted for different DM mass. For large m S , we can see that the constraints on λ HS are proportional to m χ . The reason is that the one-loop induced Higgs coupling is roughly proportional to y 2 λ HS m /m 2 S for large m S expansion. Since the relic abundance fix the combination y 2 m χ /m 2 S to be constant, the sensitivity for λ HS from Higgs precision measurement is proportional to m χ . Different from Higgs invisible branching ratio, there is no phase space suppression for m χ ∼ m h /2. One can see that for m χ 20 GeV, the Higgs-tau coupling precision measurement is the most sensitive among the three panels in the figure, while for the intermediate mass m χ between 20 to ∼ 50 GeV the Higgs invisible branching ratio measurement is better. For m χ close to m h /2, the Higgs-tau coupling measurement becomes better again due to the phase space suppression in the Higgs invisible decay.
In the middle panel of Fig. 15, we show the limits from δκ µ < 8.7%. The results from muon coupling measurements are fully analogous to tau coupling. The sensitivity is worse by a constant factor from δκ τ /δκ µ , reflecting the fact that more taus are produced due to larger Higgs-tau coupling.
Finally, it should be mentioned that the limits on λ HS from Br(h → γγ) and σ(Zh) are also very powerful as shown in section 3.3 and are able to exclude most of the parameter space for GW detection [100]. Such constraints are independent of the DM Yukawa coupling y and therefore, are complementary with the limits from h invisible and leptonic decays.

Conclusions
The GW detection opens a new window to the FOPT and the Higgs precision measurement is an inevitable path after the Higgs discovery. In this paper, we study their interplay in a specific DM model, namely lepton portal DM model. We emphasize the Higgs portal coupling in this model, which is neglected in the previous literature. The impact is investigated in two aspects. In the cosmological aspect, we have studied the parameter space allowing a FOPT and yielding detectable GW signals at the future detectors, taking LISA as an example. In the particle aspect, we have considered various new channels to further test this model: • pp → S + S − at the LHC, which mainly probes m S and m χ since the production is dominated by the Drell-Yan process, can also test the Higgs portal coupling between h and S ± (i.e. λ HS ) through the gluon-gluon fusion process.
• e + e − → S ± S ∓( * ) at the future lepton colliders, which can fill in the gaps between the LEP and LHC constrains (100 GeV m S 150 GeV), and probe the y coupling via the off-shell production of S ± .
• Exotic decays of h and Z at the future lepton colliders, which probe the couplings λ HS and/or y for the low m χ region.
• Higgs precision measurements for invisible decay branching ratios and leptonic coupling originated from one-loop contributions, which can provide the best sensitivity for the combination y 2 λ HS or λ HS assuming y satisfies the DM relic abundance requirement.
• Electron g − 2 experiments have recently came up with two sets of results. For ∆ã e , DM mass should be larger than 1 GeV, while for ∆a e , DM mass between 0.2 ∼ 2 GeV is preferred.
In summary, the future Higgs precision measurements can effectively interplay with GW detection, since they both rely on the Higgs portal coupling. The Higgs portal is allowed by this model and can contribute to the Higgs couplings to DM and SM leptons at one-loop level. Therefore, most of the GW detectable parameter space can be crosschecked by the Higgs precision measurement. It shows the rigorous interplay between the future Higgs precision measurement program and the GW detection program. Specific to the lepton portal DM model, which is hard to probe through DM direct and indirect detections, the Drell-Yan production of charged scalar pair is the useful way to probe this model but only constrains the mass parameter of the scalar and DM. We studied the Higgs mediated S ± pair production, exotic decays of h/Z and electron g − 2 experiment, which can help extending the constraints on mass parameters and also providing useful constraints on the Yukawa and scalar portal couplings.
of China under Grant No. 12005009. KPX is supported by the Grant Korea NRF-2019R1C1C1010050. KPX would like to thank the hospitality of the University of Chicago where part of this work was performed.