Leptophilic dark matter from gauged lepton number: phenomenology and gravitational wave signatures

New gauge symmetries often appear in theories beyond the Standard Model. Here we study a model where lepton number is promoted to a gauge symmetry. Anomaly cancellation requires the introduction of additional leptons, the lightest of which is a natural leptophilic dark matter candidate. We perform a comprehensive study of both collider and dark matter phenomenology. Furthermore we find that the model exhibits a first order lepton number breaking phase transition in large regions of parameter space. The corresponding gravitational wave signal is computed, and its detectability at LISA and other future GW detectors assessed. Finally we comment on the complementarity of dark matter, collider and gravitational wave observables, and on the potential reach of future colliders.


JHEP02(2019)048
including three generations of right-handed neutrinos carry unit charge, whereas the other SM fields are neutral. Lepton number is spontaneously broken by an SM singlet scalar field, giving mass to the lepton number gauge boson. Additional fermionic fields are added to cancel gauge anomalies. These additional fields are vector-like under the SM gauge group.

Gauge sector
The model is based on the gauge group SU(3) c ⊗ SU(2) W ⊗ U(1) Y ⊗ U(1) . Omitting QCD, the gauge sector of the Lagrangian is given by where W a andB are the gauge bosons of the SM weak and hypercharge gauge group respectively, andẐ is the lepton number gauge boson. The 2B µνẐ µν term leads to kinetic mixing between the hypercharge and lepton number U(1) gauge bosons. The kinetic terms can be diagonalized by a GL(2, R) transformation [40]  The hats denote fields in the gauge basis and the unhatted fields are in the basis where the kinetic terms are diagonal and canonically normalized. The model further features two scalar fields: the SM Higgs doublet transforming under SU(2) W ⊗ U(1) Y ⊗ U(1) as H ∼ (2, 1/2, 0), and a complex scalar Φ ∼ (1, 0, L Φ ) which is an SM singlet with lepton number L Φ . We let both fields acquire a vacuum expectation value (VEV) given by H = (0, v H / √ 2) and Φ = v Φ / √ 2, thus breaking the electroweak and lepton number gauge group SU(2) W ⊗ U(1) Y ⊗ U(1) to U(1) EM electromagnetism. The gauge bosons then obtain masses from the kinetic terms of the scalar fields, with the covariant derivative given by Here, g 2 , g 1 and g are the gauge couplings of the SU(2) W , U(1) Y and U(1) gauge groups respectively. The W mass is the same as in the SM, m W = 1 2 g 2 v H , whereas the mass matrix for the remaining gauge fields in the kinetic eigenbasis (W 3 , B, Z ) is given by The upper-left 2 × 2 submatrix is diagonalized by the SM weak mixing angle. If = 0, the resulting Z SM boson is mixing with the Z boson. The kinetic eigenstates are related to the physical mass eigenstates by    the would-be Nambu-Goldstone bosons G ± ,Ĝ 0 andω 0 become the longitudinal degrees of freedom of the W ± , Z and Z gauge bosons. The mass matrix for the remaining scalars is The Higgs portal term λ p H † H Φ † Φ induces a mixing between theĥ andφ fields. The mass eigenstates are defined by with the corresponding masses where we eliminated µ 2 H and µ 2 Φ using the condition that the potential (2.9) has a minimum forĥ = v H andφ = v Φ . Here, h is the SM-like Higgs with m h = 125 GeV, and φ is the lepton number Higgs which will typically have a mass m φ > m h due to the VEV hierarchy imposed by LEP constraints (see section 4.1).
In order to write Yukawa terms for the additional fermions involving the lepton number breaking scalar Φ, L Φ = 3 must be chosen. The Yukawa sector is then given by . (2.15) Note that specific values of L allow for additional terms in the Lagrangian. For example if (L , L ) = (1,4) or (−2, 1), the dark fermions can mix with the SM ones, potentially leading to flavor changing neutral currents and threatening the stability of our DM candidate. Table 1 lists the lepton number charges that allow additional terms. We will exclude these choices in the following. For any other pair of real numbers with L = L + 3, the Yukawa interactions of the exotic leptons are fully described by (2.15).
After spontaneous symmetry breaking, mass terms for the additional fermions are generated, (2.16) with the mass matrices given by (2.17) The matrices can be diagonalized via singular value decomposition, yielding the diagonal matrices M ν D = U ν † L M ν LR U ν R and M e D = U e † L M e LR U e R , where U a C are unitary matrices. For simplicity, and to avoid CP violating phases, let us assume that the Yukawa couplings c i and y  fermions combine to two charged (e 4 and e 5 ) and two neutral (ν 4 and ν 5 ) Dirac fields, which are given in terms of the original fields by cos α e sin α e − sin α e cos α e E L e L + cos β e sin β e − sin β e cos β e E R e R . (2.18) The right-and left-handed fields mix with different mixing angles unless we choose y ν = y ν and y e = y e . Thus, the resulting fermions in general are chiral with respect to both the SM and U (1) .
In the absence of the Yukawa terms (2.15), the Lagrangian exhibits a global [U(1)] 6 symmetry at the classical level, consisting of a U(1) symmetry for each additional lepton field in (2.14). The Φ Yukawa terms break this to three U(1) symmetries (one for the doublets, one for the charged singlets, and one for the neutral singlets), whereas the H Yukawa terms (in the absence of the Φ Yukawas) break the symmetry to U(1) L ⊗ U(1) L . Hence, small values of c i 1 and y ( ) i 1 are technically natural, rendering vector-like masses c i v Φ v Φ . Similarly, y SM ν i 1 are technically natural. As we will see later, v Φ 2 TeV, therefore we typically have c i v Φ y ( ) i v H . Consequently, the mixing angles α e/ν and β e/ν are usually small, and the masses are approximately given by m To simplify the discussion of the model we will restrict to the case of symmetric mass matrices, i.e. y e/ν ≡ y e/ν = y e/ν , in the following, so that α e/ν = β e/ν . Let us further assume that c e = c . The masses are then given by

JHEP02(2019)048
In particular, m ν 4 = (m e 4 + m e 5 )/2 if y ν v H c l/ν v Φ , and e 4 and e 5 are maximally mixed with α e = β e = π/4. Note that the SM neutrino masses in our model are pure Dirac masses generated from small Yukawa couplings to the SM Higgs doublet. Majorana mass terms are forbidden by the lepton number gauge symmetry, whereas mass terms arising from mixing with the exotic leptons would spoil the DM stability and are therefore avoided by suitable choices of the lepton number charges.

RG running
Before exploring the phenomenology of the model, let us first consider the renormalization group running of the lepton number gauge coupling g .
The running of a gauge coupling g is governed by the beta function For a U(1) gauge group, the one-loop beta function is where the sums run over all Weyl fermions and complex scalars charged under the gauge group with charge Q f /s , respectively. For the lepton number gauge group we get contributions from the SM leptons (with unit charge), the two additional generations of vector-like leptons (with charge L and L ), and the lepton-number-breaking scalar (with charge L Φ ). Thus, where we used the lepton number charges L = L + 3, L Φ = 3 and the number of SM flavors N f = 3. The dependence of the gauge coupling on the scale µ is consequently given by where g 0 = g (µ 0 ) and b = 35 + 16L + 16 3 L 2 . U(1) has a Landau pole when the second term in the bracket in equation (2.24) is of order unity, i.e. at the scale µ = Λ with Figure 1 shows the Landau pole Λ normalized to the scalar VEV v Φ as a function of g 0 = g (m Z ) for different values of the charge L . Certainly, we want the Landau pole to occur significantly above the Z mass and above v Φ , otherwise the validity of our perturbative results would be questionable. This requires choosing g (m Z ) 0.5 for most values of L . The slowest running is obtained for L = −3/2, which we however excluded since it allows for Majorana mass terms.
To prevent the gauge coupling to run into a Landau pole at low scales, we choose L = −1/2 in the remainder of this paper. 2 In this case g (m Z ) ≈ 1 is acceptable, with the Landau pole located almost two orders of magnitude above v Φ . For v Φ = 2 TeV this implies an upper bound of m Z 6 TeV for the mass of the Z boson. For most of the considerations that follow, the exact value of L merely matters anyway. An exception are the dark matter observables considered in section 3, where we therefore also consider different values for L .

Leptophilic dark matter
Provided that we avoid the specific choices of L and L discussed in section 2.3, the model features a global U(1) L +L symmetry under which all SM fermions are neutral whereas the additional leptons have unit charge, and which is free of anomalies. This symmetry persists when the EW and U(1) gauge symmetries are broken and ensures the stability of the lightest dark lepton. If neutral, it is a candidate for dark matter.
For the remainder we identify ν 5 as the DM candidate (which can always be achieved by defining the mixing angles in (2.18) accordingly) and relabel it as ν DM ≡ ν 5 . Since direct detection experiments exclude DM candidates with unsuppressed couplings to the SM Z boson, ν DM should be composed predominantly of the SM singlets ν L and ν R . Consequently α ν and β ν should be small. In particular, this requires c ν < c .  Table 2. Default values for the model parameters (assuming y ν/e = y ν/e and c e = c l ) used throughout this paper, unless specified otherwise. For negligible sin(θ DM ), the mostly-doublet, heavy neutrino mass is given by m ν4 (m e4 + m e5 )/2 = 1.75 TeV. at least one of y ν and y ν should be non-vanishing, otherwise an additional global U(1) symmetry remains unbroken and the next-to-lightest dark sector state (either ν 4 or e 4/5 ) would be stable as well.
The model is implemented in FeynRules [41] and subsequently mircOMEGAs [42] was used to calculate the relic density as well as direct and indirect detection constraints. Unless specified otherwise, we use the parameters listed in table 2. We here relabelled θ DM ≡ α ν = β ν .

Relic abundance
Assuming that ν DM is a thermal relic, its abundance is predominantly set by its annihilation cross section to two leptons through a s-channel Z (figure 2(a)). Other possible annihilation channels are annihilation to gauge or scalar bosons through an intermediate h or φ, or to fermions via a Z boson. The former is suppressed by the Higgs-scalar mixing, whereas the latter can arise from Z − Z mixing or by an admixture of the SM component in the DM. The doublet-singlet mixing further allows for t-channel annihilation to two bosons, and a small mass splitting between the DM and the other exotic leptons can lead to coannihilation. See [4] for more details.
The parameter regions in the m DM − v Φ and m DM − m Z planes that reproduce the DM relic abundance of h 2 Ω DM = 0.1198 ± 0.0015 measured by the Planck satellite [43] are shown in figure 3 for different values of L . We assume a lepton number gauge coupling of g = 0.1 and a scalar self-coupling of λ Φ = 0.5 as well as Yukawa couplings c = 1.5 and y e = 0 when scanning over the VEV (left panel), and a scalar VEV of v Φ = 2 TeV when varying the Z mass (right panel). The remaining parameters are set to the values specified in table 2.  The colored regions yield a DM abundance that lies within two standard deviations around the Planck measurement. For each value of m Z we typically obtain two viable values for the DM mass, one below and one above the m DM = m Z /2 resonance. In figure 3(a), we can in addition also see the scalar resonance at m DM = m φ /2. To guide the eye, the resonances are indicated by dashed, dark gray lines.
Since the Z predominantly decays into ν DM and the other heavy leptons, its width increases with L . For L = −1/2, the Z is rather narrow and the DM mass is restricted to values close to half of the Z mass. For larger charges, the resonance becomes broader and the DM mass can be lower. The light gray, dashed line in figure 3(b) indicates the regions in which the width of the Z exceeds 10 % of its mass for L = 3/2.
The dependence of the DM relic density on the masses of the non-DM heavy leptons e 4 , e 5 , and ν 4 is shown in figure 4, assuming that they all have the same mass m HL . The colored regions again yield the measured DM abundance, now assuming L = −1/2. The colors correspond to a relative mass splitting ∆ m ≡ (m HL − m DM )/m DM of 1 % (blue), 2 % (red), 5 % (green), and 10 % (purple) between the DM and the heavy leptons.
For high DM masses, the heavy lepton masses affect the relic density only by changing the Z width. However, for lower DM masses the abundance is no longer set by annihilation of the DM to SM leptons as depicted in figure 2(a), but via co-annihilation. In this case, the heavy lepton (HL) abundance is depleted by annihilation of e 4 , e 5 and ν 4 to SM particles through electroweak processes as shown in figure 2(b). This depletion is transferred to the DM abundance by the electroweak processes depicted in figure 2(c) in which a DM particle scatters off SM particles and changes into a heavy lepton. These processes require sin θ DM = 0, but we need this assumption anyways to ensure that there is only a single DM component. As the diagram 2(b) is dominanted by electroweak processes, the relic density in this regime is independent of the Z mass. Figure 4 assumes sin θ DM = 0, however modifying the mixing within the range allowed by direct detection (see section 3.2) does not alter the result.
Varying the remaining parameters of the model only has a minor effect on these results. The scalar mass m φ and the h − φ mixing angle θ H only have an effect in the region of the φ resonance m DM = m φ /2 (or the h resonance).

Direct and indirect detection
Direct detection experiments strongly constrain DM couplings to the SM Z boson via scattering off nuclei. For small values of the kinetic mixing parameter , the coupling is given by where s W = sin θ W , s ξ = sin ξ, s DM = sin θ DM and c DM = cos θ DM . The first term originates from the heavy doublet neutrino N = N L + N R (which has vector-like couplings to the SM Z) mixing into the DM, whereas the second part comes from the chiral DM − Z coupling. The axial part of the latter is also modified by the DM mixing via the ν 4 ν 4 Z vertex, the vector part remains untouched by θ DM since it here enters as (c 2 DM + s 2 DM ). Figure 5 shows the constraints on the Z − Z and DM mixing as a function of the Z mass, obtained from direct detection limits on spin-independent DM-nucleus scattering. The solid lines correspond to current constraints from the XENON1T experiment based on  one tonne times year of data acquisition [44], the long dashed lines indicate the prospective sensitivity of LZ [45], and the dash-dotted lines show the projected reach of DARWIN [46]. At each parameter point the DM mass is fixed to a value reproducing h 2 Ω DM = 0.1198 (chosing the value below the Z resonance), the charges are taken to be L = −1/2 (blue) or L = 3/2 (red), and the VEV is set to v Φ = 2 TeV.
Current direct detection experiments can probe kinetic mixing parameters in the percent range and DM mixing angles of sin θ DM ∼ 0.015 − 0.025, depending on the Z mass. With LZ, DM mixing angles of sin θ DM ∼ 0.006−0.01 and kinetic mixing in the sub-percent range can be reached. DARWIN can prospectively exclude sin θ DM ∼ 0.004 − 0.006, and sub-per-mill kinetic mixing for m Z 1 TeV. The constraints for L = 3/2 are stronger than for L = −1/2 since the latter case leads to higher DM masses, whereas the former case gives DM masses below 500 GeV.
Beside Z-mediated DM-quark interactions, nuclear scattering can also proceed via Higgs (or φ) exchange, either through Higgs-scalar mixing or by direct DM-Higgs couplings. However, since the Higgs only weakly couples to nuclei, the direct detection constraints on the Higgs mixing angle are much weaker than on θ DM or ξ.  lower φ masses on the other hand, the scattering cross section is reduced due to φ − h interference effects, leading to weaker direct detection limits. The corresponding constraints are shown in figure 6.
A further, indirect way of probing DM is via its annihilation to SM particles. Since the observation of charged particles suffers from large uncertainties associated with their propagation through the Galactic halo, we here only consider indirect detection constraints from γ-ray searches. As photons travel unperturbed by Galactic magnetic fields, they can be traced back to their production site, allowing for constraints on DM annihilation by observing photons from regions with a high DM density.
In our model the DM typically annihilates through the lepton number gauge boson Z into SM leptons with equal branching ratios. Photons are thus predominantly produced as secondary products from annihilation to charged leptons. Direct production of monochromatic photons is possible via annihilation through the scalar bosons h and φ, however, this is suppressed unless resonant.
We tested the annihilation of thermally produced DM (i.e. satisfying the relic density constraint) in our model against current limits from observations of dwarf spheroidal galaxies by MAGIC and Fermi-LAT [48], and of the inner Galactic halo by H.E.S.S. [49], as well as against γ line searches from Fermi-LAT [50] and H.E.S.S. [51]. The strongest limits come from secondary produced photons from annihilation into τ leptons. However, the current sensitivity reaches the level of the annihilation cross section required for thermal production (which in addition is reduced by the branching ratio of 1/6 into tauons) -13 -

JHEP02(2019)048
only for DM below 100 GeV, which, even for light Z masses, is below the dark matter masses predicted by our model (cf. figure 3). This also holds for the projected sensitivity of Fermi-LAT, assuming a 15-year data set of 60 dwarf spheroidal galaxies [52]. On the other hand, a next-generation γ-ray observatory such as the Cherenkov Telescope Array (CTA) [47] will be able to exclude Z masses between roughly 670 GeV and 1.46 TeV if L = 3/2. The corresponding limit is indicated by the vertical, dashed lines in figure 5 In principle, our model can furthermore be probed through its neutrino sector. Modifications of the neutrino interactions with the SM leptons arise from Z exchange or kinetic mixing, and DM-neutrino interactions can be mediated by a Z or Z boson. Neutrino couplings to the lepton number breaking scalar and SM Higgs boson are suppressed by the neutrino Yukawa couplings. For the range of Z and DM masses considered here, no constraints are obtained from current data [53,54].

Collider phenomenology
While new physics that couples directly to quarks and gluons is nowadays severly constrained by direct searches at the LHC, the situation is different for the leptophilic new physics model we are considering here. In this model constraints predominantly arise from a combination of LEP limits as well as direct and indirect LHC measurements, such as e.g. Higgs data. In the following we present an overview of the most important constraints on the lepton number gauge boson, the extended Higgs sector and the new leptons introduced in our model, and comment on the prospects for detection at the high-luminosity LHC and future colliders.

Z constraints
As in the absence of kinetic mixing the lepton number gauge boson does not couple to quarks, the strongest constraints on the Z boson come from LEP II. These exclude Z masses below the maximal LEP center-of-mass energy of 209 GeV (except for tiny gauge couplings g < 10 −2 [55,56]) and severely restrict the lepton number breaking VEV v Φ through 4-lepton contact interactions.
A heavy Z induces effective contact interactions between electrons and charged (SM) leptons. Since the Z couplings to SM fermions are vector-like, the corresponding contact interactions are given by (neglecting kinetic mixing) which interfere destructively with the SM Z boson for center-of-mass energies above the Z-pole. LEP puts a 95 % lower bound of Λ > 20 TeV [57] on the scale suppressing this contact interaction, related to the model parameters by Λ 2 = 4π To be conservative we set the VEV to v Φ = 2 TeV in the following.  Future e + e − colliders have the potential to substantially tighten these bounds. For instance, the ILC with a center-of-mass energy of √ s = 1 TeV can basically exclude Z masses below 1 TeV (unless the gauge coupling is below g 7.6 × 10 −8 ), and constrain the VEV to be above v Φ 15 TeV at 90 % CL using muon contact interactions [17].
At the LHC, the Z is rather hard to produce. In particular, in the absence of kinetic mixing it is predominantly produced by pair-producing SM leptons that radiate off a Z . The Z can then be detected from its decays to charged SM leptons.
To obtain a rough estimate of the detection prospects, we calculate the parton level cross section for pp −→ Z with CalcHEP 3.6 [58], where can be any SM lepton, including neutrinos. For the decay we use the narrow width approximation, assuming that the Z decays into SM leptons only. In this case, the corresponding branching ratio to charged leptons is Br (Z → + − ) = 50 %. The cross section as a function of the Z mass is shown in figure 7 for the LHC with center-of-mass energies of 13 TeV (blue) and 14 TeV (red), as well as for a 100 TeV collider (green). The gray lines indicate the cross sections that would produce 10 Z s assuming integrated luminosities of 300 fb −1 and 3 ab −1 .
With the current LHC data, no constraints can be put on the Z mass if kinetic mixing is absent. With 300 fb −1 , the LHC would not produce a sufficient amount of Z s. At the high luminosity LHC, Z masses below ∼ 400 GeV can be reached. The prospects for a 100 TeV collider are more promising. With 3 ab −1 , more than 10 events are produced for masses up to 2.5 TeV, extending the reach into the multi-TeV region.
In the presence of kinetic mixing the situation is different. The lepton number gauge boson then couples to quarks with couplings proportional to the kinetic mixing parameter and the quark hypercharge. It can thus be produced directly in proton-proton collisions and be searched for as a dilepton resonance, giving constraints in the m Z vs. plane. (dotted purple) are also shown. Currently, kinetic mixing can be probed in the percent range, the HL-LHC can prospectively reach the sub-percent range for light Z . Again, the cross sections have been calculated with CalcHEP [58], assuming a narrow Z width with a branching ratio of 1/3 to light leptons (e or µ). The projections have been obtained assuming that the limits on cross section ratios provided by CMS [60] do not change when increasing the center-of-mass energy from 13 TeV to 14 TeV, and that the exclusion reach scales with the square root of the luminosity.
The kinetic mixing can also be probed via its effects on SM precision measurements at electron-positron colliders [61]. However, as these effects are suppressed for high Z masses, the LHC provides the strongest constraints in the mass range considered here.

Higgs constraints
The scalar sector of our model is subject to constraints from measurements of the properties of the 125 GeV Higgs boson at ATLAS and CMS, as well as from null-results of searches for scalar bosons at different masses.
In our model, the SM Higgs properties can be modified by three effects: the h−φ mixing which modifies all SM Higgs couplings, modifications of Higgs couplings to electroweak gauge bosons (in particular to two photons) by loops of heavy charged leptons, and decays to BSM states (if kinematically accessible). However, given the lower bound on the leptonnumber-breaking VEV (4.2), the new states are typically too heavy for the SM Higgs to decay into, so that the last effect is absent in most of the parameter space.
The mixing between the lepton-number-breaking scalar and the SM Higgs boson given by equation ( Neglecting additional Higgs decay channels and further modifications of loop-induced Higgs couplings discussed below, the production cross sections are modified by a factor cos 2 θ H , whereas the branching ratios remain unchanged as the cosine factors in the partial and total widths cancel. Thus, the signal strengths are µ = cos 2 θ H . An estimate of the limit on the mixing angle can be obtained from the global signal strength. The current CMS measurement is µ = 1.17 ± 0.10 [62]. This gives a 95 % exclusion of Loops of the dark electrons e 4 and e 5 can contribute sizeably to the h −→ γγ and h −→ Zγ rate. In the SM, these rates are given by [63] where α is the electromagnetic coupling constant, The expressions for the form factors A s for a spin s particle running in the loop can be found in [63]. The sums run over all charged fermions that couple to the Higgs. N f c is the colorrepresentation of the fermion, Q f is its electric charge, andv f is the fermion's (reduced) vector-coupling to the Z boson. In the SM, the dominant contribution from fermions comes from the top quark with N t c = 3, Q t = 2/3, andv t = 1 − 8 3 s 2 W . Equations (4.5) and (4.6) assume that the fermions couple to the Higgs with a vertex factor proportional to their masses. The top quark in the SM for instance has a vertex factor −i yt for the SM-like Higgs. 3 The correct result can then be obtained by rescaling the heavy lepton contributions by a factor , where s = h, φ and f = e 4 , e 5 . Due to the scalar 3 The corresponding interactions of φ are obtained by replacing cH → sH and sH → −cH . mixing, the SM contributions further get a factor of c H or s H for h and φ, respectively. We thus obtain . (4.9) The corresponding width for the φ scalar can be obtained by replacing The leading QCD corrections can be included by multiplying the top contribution by 1 − αs π [63]. We evaluated constraints from direct Higgs searches using If the new leptons have sizeable couplings to the Higgs boson, the Higgs signal strengths in different channels can vary due to the loop contributions to h → γγ and h → Zγ decays. Figure 10 shows the current LHC limits from [62,[72][73][74] as a function of the mixing angle and the heavy electron Yukawa couplings c and y e . If the Φ Yukawa coupling c is small, the dark electrons gain their mass predominantly from electroweak symmetry breaking and hence strongly contribute to the h −→ γγ rate. Thus, the Yukawa coupling to the Higgs doublet y e is also restricted to be small. For large c , the heavy electron contributions in (4.8) are mass-suppressed, so that y e can take large values without modifying Γ (h −→ γγ) beyond the experimentally allowed limits.
Note that for c = 0 the charged, heavy lepton masses are given by m e 4/5 = yev H √ 2 . The LEP limit on the mass (see section 4.3) then constrains the Yukawa coupling to y e > 0.57, hence the entire region allowed by h → γγ and h → Zγ is excluded in this case. Similarly, (2.20) implies |y e | < 0.24 for c = 0.1 and v Φ = 2 TeV. Further note that the exclusion for c = 10 is shown despite severly challenging the bounds of perturbativity to illustrate the constraints in the limit of large c .

Constraints on heavy leptons
Whereas the DM candidate is mostly an SM singlet, the remaining heavy leptons carry electroweak charge and can hence be produced at colliders. Direct searches for charged heavy leptons at LEP set a lower limit of roughly 100 GeV on the e 4 and e 5 mass [75]. LHC limits on the dark leptons can be obtained by recasting supersymmetry (SUSY) searches for electroweakly produced charginos and neutralinos.   Figure 11 shows the cross section for the production of two charged heavy leptons (dashed, red), heavy positrons with a heavy neutrino (solid blue), heavy electrons with a heavy anti-neutrino (dash-dotted, green), and pairs of heavy neutrinos (dotted, purple) in proton-proton collisions at √ s = 13 TeV calculated using CalcHep [58]. We here take the most optimistic scenario for heavy lepton production in which all heavy leptons have the same mass. As the DM mixing angle θ DM is restricted to be small by the direct detection constraints in section 3.2, and therefore only has a negligible effect on the production cross section, it can be set to zero.
Due to the U(1) symmetry that stabilizes the DM, the heavy leptons can only decay amongst themselves or to dark matter. Consequently, to allow the lighter dark electron to decay, θ DM = 0 is required. Otherwise the model would have a charged DM population and therefore be excluded. However, even a (negligibly) small amount of DM mixing is sufficient to let the exotic particles decay fast enough to avoid this problem.    Figure 14. CMS 95% exclusion limits from [77], assuming equal masses for all heavy leptons and that they decay within the detector.
The charged heavy leptons typically decay into DM and a (potentially off-shell) W boson. Depending on the masses, decays to other exotic leptons can also be possible, these are however suppressed by phase space. The heavy neutrino can decay to DM and a(n off-shell) Z boson or, if m ν 4 > m DM + m h , to DM and an SM Higgs boson. In the latter case, the branching rations to DM + Z and DM + h are roughly 50 %. Figure 12 shows the branching ratios of ν 4 as a function of the mass for a DM mass of 100 GeV and a mixing angle of sin θ DM = 0.01.
At the LHC, the heavy leptons can be searched for by looking for missing transverse energy in association with W , Z, or h bosons (or SM lepton pairs in the off-shell case). These searches have been performed by ATLAS [76] and CMS [77] in the context of simplified SUSY models, using 36 fb −1 of data recorded at the LHC with √ s = 13 TeV.
They assume that the lightest neutralinoχ 0 1 is the lightest supersymmetric particle (LSP) and consider the process pp −→χ ± 1χ 0 2 , where the lightest charginoχ ± 1 decays to the LSP plus a W boson, and the next-to-lightest neutralinoχ 0 2 to the LSP plus Z or h. The respective 95 % exclusion bounds from CMS can be found in figure 8 of [77].

JHEP02(2019)048
The production cross section for the corresponding process in our model is depicted in figure 13, again assuming m e 4 = m e 5 = m ν 4 (sin θ DM = 0). The respective exclusion bounds from CMS [77] are shown in figure 14, taking the limits assuming Br χ 0 2 →χ 0 1 Z = 100% for m DM < m e 4 /e 5 /ν 4 < m DM + m h , and Br χ 0 2 →χ 0 1 Z = Br χ 0 2 →χ 0 1 h = 50% for m e 4 /e 5 /ν 4 > m DM + m h . The LHC can currently exclude heavy lepton masses below m e 4 /e 5 /ν 4 200 GeV and DM masses below m DM 150 GeV. For the co-annihilation region discussed in section 3.1, the mass splitting between the charged states and the DM candidate becomes very small, so that the searches used here become inefficient. Instead it has been shown that these regions can be probed by mono-jet, mono-Z and disappearing track searches, with masses of up to 200 GeV reachable at the LHC and up to 1 TeV at a future hadron collider [78][79][80][81][82][83].

The lepton number breaking phase transition
In the early Universe, spontaneously broken symmetries are typically restored due to thermal effects. This can be seen by computing finite-temperature corrections to the effective potential of the scalar fields whose vacuum expectation values break the symmetries. In the model at hand, these are the scalarφ breaking the U(1) lepton number gauge symmetry, and the Higgs fieldĥ which breaks SU(2) W ⊗ U(1) Y . As a consequence, a symmetry breaking phase transition (PT) occurred in the history of the Universe.
At high temperatures, the global minimum of the finite-temperature effective potential is at the origin, i.e. SU(2) W ⊗ U(1) Y ⊗ U(1) is unbroken. When the temperature drops, the potential changes and at some point develops a minimum at non-vanishing field values. Whether both fields develop non-zero VEVs at the same time or independently at different temperatures depends on the parameters of the model. Since the portal interaction between the two scalars is restricted to be small (see section 3.2 and 4.2), and due to the hierarchy of the VEVs (v H = 246 GeV, v Φ 1.9 TeV), the lepton number breaking PT typically occurs first at temperatures in the TeV range, leaving the EW symmetry unbroken. EW symmetry breaking subsequently occurs as in the SM at a temperature of T 160 GeV [84].
The phase transition between the unbroken and the broken phase can proceed in basically two different ways. In a first-order phase transition the scalar fields tunnel or fluctuate through a potential barrier between the two phases, whereas in a cross-over transition they smoothly evolve from one vacuum to the other. The EWPT in the SM falls into the latter category [84].
Here we are mainly interested in transitions of the first kind. These have a much richer phenomenology, potentially giving rise to a stochastic background of gravitational waves as discussed in section 6, or providing out-of-equilibrium dynamics which are, along with baryon number, C and CP violation, required by the Sakharov conditions for baryogenesis [85].
A cosmological first-order phase transition can occur when at some critical temperature T c the scalar potential develops two degenerate minima separated by a barrier. The transition then proceeds through the nucleation of bubbles of the true vacuum (i.e. the global minimum at low temperatures) in the sea of the false vacuum (i.e. the global mini--22 -JHEP02(2019)048 mum at high temperatures). The energy released from the potential difference of the two vacua drives the expansion of the bubbles, which collide and merge, and eventually fill the whole universe.
The bubble nucleation rate per unit volume at finite temperature is given by [86][87][88] where A(T ) ∼ T 4 is a prefactor and S E is the Euclidean action evaluated for the O(3) symmetric bounce solution. In the early Universe the process of bubble nucleation has to compete against the expansion of the Universe, given by the Hubble rate H. Furthermore, the surface energy would make nucleated bubbles collapse unless the energy gain from the potential difference is larger. Consequently, the transition does not occur immediately at the critical temperature, but at a lower temperaturethe so-called nucleation temperature T n . It is defined by the temperature at which the probability to nucleate one sufficiently large bubble per Hubble volume is of order one.
The bubbles nucleated within one Hubble patch proceed to expand and collide, until the entire volume is filled with the true vacuum.

Finite temperature effective potential
The daisy-resummed one-loop effective potential of the scalar fields at finite temperature is given by [87,89,90] Here, h and φ are classical background fields, 4 and V 0 is the tree-level potential of h and φ given by V CW is the Coleman-Weinberg potential [91] which includes all temperatureindependent one-loop corrections at vanishing external momentum. In dimensional regularization with MS renormalization it reads The sum runs over all particles that couple to the scalars. Since we compute the potential in Landau gauge this also includes Goldstone bosons. Here n i are the numbers of degrees of freedom of each particle, which include an additional factor −1 for fermions, and m 2 i (h, φ) are the field dependent masses. The constants C i are C i = 3 2 for scalars or fermions and C i = 5 6 for gauge bosons.

JHEP02(2019)048
V c.t. includes the non-MS part of the counter-terms, i.e. the difference of the counterterms to the usual MS counter-terms. These can be derived from the Coleman-Weinberg (and tree-level) potential and take the form The exact form of the counter-term couplings depends on the renormalization conditions imposed.
The temperature-dependent one-loop corrections to the potential are contained in where the thermal functions for bosons (J B , "−" case) and fermions (J F , "+" case) are For T 2 m 2 , these can be expanded as [87,89,90] where a b = 16π 2 exp(3/2 − 2γ E ) and a f = π 2 exp(3/2 − 2γ E ). Finally, the resummed daisy corrections are wheren i indicates that the sum only runs over scalars and the longitudinal gauge boson degrees of freedom, and Π i (T ) are the leading thermal mass corrections. In figure 15 we show the effective potential at different temperatures, to illustrate the individual steps of the symmetry breaking process. The model parameters are given in the figure caption. Counter-terms are set by the condition that the first and second derivatives of the zero-temperature effective potential are the same as for the tree-level potential.
At high temperatures, the global minimum of the effective potential is in the symmetric (unbroken) vacuum (h, φ) = (0, 0). As the Universe cools down, a second minimum starts to form at non-vanishing values of φ. At T c 835 GeV, the two minima are degenerate. At lower temperatures, the second minimum (h, φ) ∼ (0, 1.1 TeV) is the global minimum and breaks lepton number, whereas the electroweak symmetry remains unbroken. This minimum is separated from the symmetric minimum by a potential barrier. Thus, to transition to the global minimum, the field has to tunnel (or remain in the symmetric vacuum until the barrier disappears).
As the Universe cools further, the minimum at the origin disappears and the global minimum moves towards the zero-temperature lepton number breaking VEV (h, φ) = (0, 2 TeV). Subsequently, at T 160 GeV, the minimum starts to shift to non-vanishing Higgs field values, breaking the electroweak symmetry in a cross-over transition. Eventually, the Universe ends up in today's vacuum (h, φ) (246 GeV, 2 TeV).

A first-order lepton number breaking PT
In this section, we examine the lepton number breaking phase transition in the limit of negligible portal coupling λ p between the SM Higgs and the scalar Φ. Further assuming that the kinetic mixing of the gauge bosons as well as the exotic Yukawa couplings c i -25 -

JHEP02(2019)048
and y i of the heavy leptons are small, we can study a simplified version of the effective potential in which only the lepton number breaking scalar and the lepton number gauge boson are considered.
In this case, the tree-level potential simplifies to Setting the lepton number breaking VEV to v Φ = 2 TeV, in agreement with the LEP constraint, the model is therefore fully specified by m Z and m φ . The field dependent masses of the scalar, the gauge boson, and the Goldstone boson are given by The thermal mass corrections are (Π φ = Π ω 0 = Π Φ ) where the first part of Π Z L comes from the scalar and the second part from the SM and exotic leptons. The subscript L of Π Z L indicates that only the longitudinal part of the Z boson receives a thermal correction. We further use an on-shell scheme, imposing the conditions This ensures that the VEV and the scalar mass at zero temperature remain at their treelevel values. Here, ∆Σ ≡ Σ(m 2 φ ) − Σ(0) is the difference of the scalar self-energy evaluated at the tree-level mass and at zero-momentum, see appendix A. The second derivative of the Coleman-Weinberg potential in the vacuum suffers from logarithmic divergences originating from the vanishing Goldstone masses. These are IR divergencies and are due to the fact that the effective potential is evaluated at vanishing external momentum. However, the scalar self-energy at zero-momentum suffers from the same divergences, hence its presence in the second condition above [92]. The divergences in ∆Σ and ∂ 2 V CW /∂φ 2 then cancel, ensuring that we obtain finite counter-terms.
We use the numerical package CosmoTransitions [93] to evaluate the effective potential and to analyze the phase transition. Fixing the VEV to v Φ = 2 TeV (and setting the renormalization scale to µ R = v Φ ), we identify the region in the m φ − m Z parameter space at which a first-order phase transition occurs.
In this model, the potential barrier between the vacua is generated by thermal corrections from gauge boson loops (note the cubic term in the high-T expansion of the bosonic thermal integral (5.9)), i.e. the larger the gauge coupling (and hence also the Z mass) the higher the barrier. Increasing the scalar mass on the other hand increases the quartic coupling, which in turn reduces (the relative height of) the barrier. Thus, first-order phase transitions can be obtained for m Z m φ ; strong transitions occur for m Z 2m φ . The term "strong" here refers to transitions in which the VEV (or more precisely the distance between the two degenerate minima in field space) at the critical temperature is larger than the critical temperature itself, i.e. v Φ (T c )/T c 1. This measure is often employed in the context of baryogenesis [87]. Although the renormalization conditions (5.15) ensure that the zero-temperature potential has a minimum at φ = v Φ , this minimum is not necessarily the global minimum. In particular, if the gauge boson mass m Z is much bigger than the scalar mass m φ , the potential develops a global zero-temperature minimum at φ = 0, i.e. the Coleman-Weinberg corrections restore the symmetry already at T = 0. 6 This is the case in the white area labeled by " φ 0 = 0" above the colored region in figure 16(a) (and above the dotted line in figure 16(b)), which is of course excluded since it would imply the existence of a second 5 Note that isolated, white dots in the colored regions do not indicate that particular parameter points do not feature a first-order PT, but are due to numerical artifacts. 6 Of course, the physical scalar and Z masses become m φ = 0 and m Z = 0 in this region. Hence, the x and y axes should be interpreted as 2λΦv 2 Φ and 3g vΦ respectively, where vΦ = 2 TeV has no physical meaning.
-27 -JHEP02(2019)048 massless gauge boson with significant couplings to leptons. Furthermore, even a global minimum at φ = v Φ does not automatically ensure that today's Universe has transitioned to the true vacuum. If the barrier is very large with a small potential difference between the two vacua, which is the case close to the region in which the potential has a global zero-temperature minimum at φ = 0, the tunneling probability is too low. Therefore the field is stuck in the false vacuum and does not tunnel. This corresponds to the parameter region labeled "no tunneling" in figure 16(b). 7 On the other hand, for m Z m φ no significant barrier is induced and there is no temperature at which the potential has degenerate minima. Also, if the potential barrier separating the phases is very shallow, it might disappear before bubbles are nucleated. In both cases the transition occurs without tunneling as a cross-over 8 and no gravitational waves are generated. This happens in the areas labeled "no barrier" or "cross-over" in figure 16.
So far, to simplify the parameter space to two dimensions we neglected the contributions from the dark Yukawa couplings of the fourth and fifth generation leptons to the effective potential. However, if the leptons are heavy, the Yukawa couplings are large and the potential can be modified significantly. This is in particular the case for large Z masses, where the exotic leptons are required to be heavy in order to obtain the correct DM relic abundance.
For simplicity, we here again assume that the exotic electrons and the exotic neutrino have equal masses, m HL ≡ m e 4 = m e 5 = m ν 4 , i.e. that the SM Higgs Yukawa couplings y ν and y e in Equation (2.15) vanish. The field dependent masses are then given by The Yukawa couplings further contribute to the scalar thermal mass correction (5.14), which becomes (5.17) Figure 17 shows the values of the scalar and Z masses that lead to a first-order phase transition with the corresponding nucleation temperature, for two different choices of the DM and heavy lepton masses. The region that gives rise to a first-order transition for vanishing fermion couplings (cf. figure 16(b)) is indicated by the dashed lines.
As expected, for light dark leptons (low Yukawa couplings) the situation changes only marginally with respect to the case assuming vanishing Yukawas. However, for higher fermion masses, the region that yields a first-order PT changes, and the nucleation temperature decreases.
In the parameter region labeled " φ 0 = v Φ restored" in figure 17 (b) mDM = 500 GeV, mHL = 1 TeV Figure 17. Parameter points giving rise to a cosmological first-order phase transition with a nucleation temperature T n , including the contribution from dark matter and the exotic leptons. The dashed lines indicate the corresponding region neglecting the fermion contributions (cf. figure 16(b)).
In the gray shaded region, the potential becomes unstable below φ = 100 GeV; above the gray, dotted line it is stable up to φ = 10 6 TeV. The stars indicate benchmark points listed in appendix C.
absence of fermions. If the dark sector leptons are included, their contributions have the opposite sign and partially cancel the bosonic ones, and the global minimum at φ = v Φ is restored. Hence, the region allowing for a first-order PT is extended. On the other hand, if the fermionic corrections overcome the bosonic ones at high field values, the potential is destabilized as it is not bounded from below. This occurs for low Z and φ masses. The gray shaded regions are excluded since the potential becomes unstable at field values below φ = 100 TeV, i.e. V eff (100 TeV) < V eff ( φ 0 ) at T = 0. Above the gray, dotted curve the potential is stable even up to φ = 10 9 GeV. Note however that a reliable evaluation of the potential at such high field values requires the inclusion of renormalization group effects. At high temperatures, the loop corrections from the fermions give a positive contribution ∼ φ 2 T 2 , whereas they do not contribute to the cubic terms (note that there is no y 3 term in (5.10)). As a consequence, the finite-temperature corrections restore the symmetric minimum at lower temperatures, reducing the nucleation temperature.
Finally, to properly connect to the DM picture let us require that the DM candidate has the correct thermal abundance. 9 Figure 18 shows the nucleation temperature for the corresponding phase transitions, assuming that m HL = 1.5×m DM . At each parameter point in the m φ − m Z plane we use micrOMEGAs to find the value of the DM mass that yields the measured abundance, picking the value below the Z resonance (i.e. we are sitting on the upper branch of the blue line in figure 3). Again, the dashed lines indicate the parameter region that provides a first-order PT if the fermions are neglected.
As the DM mass required to obtain the correct abundance increases with the Z mass and is mostly independent of the scalar mass, the effects of including the vector-like leptons are stronger for larger Z masses. Hence, the fermionic corrections restore the T = 0

Gravitational waves
A first-order phase transition proceeds through thermal fluctuations or quantum tunneling by the nucleation of bubbles of the broken phase in the sea of the unbroken phase. The bubbles then expand and collide, producing a stochastic background of gravitational waves. If the phase transition is sufficiently strong first-order, these gravitational waves may be detectable by future observatories such as LISA [88,94].
In the following we will briefly review the relevant parameters and the calculation of the GW spectrum, and explain how we evaluate the detectability of the resulting spectra. See [88,[95][96][97] for some recent, more comprehensive reviews. The corresponding results are presented in section 6.4.

Transition parameters
For a phase transitions occurring at a temperature T * , the GW spectrum depends on four parameters. The most important ones are the transition temperature T * which we approximate by the nucleation temperature T n , the inverse time duration β of the transition divided by the Hubble rate 1) and the ratio of the latent heat L(T * ) released by the transition (i.e. the difference of the energy densities in the two vacua) to the radiation density where g * is the effective number of relativistic degrees of freedom.

JHEP02(2019)048
The spectrum further depends on the bubble wall velocity v w . For phase transitions occurring in a plasma (in contrast to those in a vacuum dominated epoch) the expanding bubbles are exposed to a friction force exerted by the particles of the plasma. If the friction is sufficient to compensate the driving force from the energy liberated by tunneling, the bubble walls reach a terminal velocity, otherwise they accelerate perpetually with v w approaching the speed of light. This latter case is referred to as run-away [98]. In this case, the energy transferred to the plasma saturates, and the surplus energy accelerates the scalar bubbles, giving rise to a non-negligible GW spectrum from bubble collisions (cf. section 6.2). It has however been argued that the presence of vector bosons that gain a mass in the PT produces friction from transition radiation that prevents the bubble from such a behavior [99]. Still, bubble velocities of v w ∼ 1 can be reached.
As the detectability of the generated stochastic background eventually only mildly depends on the wall velocity, we simply take v w = 1 in the following. See appendix B for a brief study of the impact of the wall velocity on the detectability.

GW spectrum
Cosmological first-order phase transitions can produce gravitational waves mainly via three mechanisms: the collision of the bubble walls themselves, acoustic production from sound shells surrounding the bubble walls, and turbulence effects in the plasma. The corresponding GW spectra can be simulated numerically or estimated analytically [100][101][102], and are then expressed in terms of quantities that can be derived from the effective potential by fitting to the corresponding results. The full spectrum is given by where Ω GW (f ) = 1 ρc dρ GW d log f is the spectral GW energy density normalized to the critical energy density ρ c = 3H 2 /(8πG), and H = h × 100 km Mpc −1 s −1 is the Hubble rate.
The first contribution h 2 Ω φ (f ) to the GW spectrum comes from the scalar bubbles nucleated during the transition. Due to the spherical symmetry of the bubbles, the potential energy that is exempted when the field tunnels cannot be released directly into gravitational waves. It can thus only drive the expansion of the bubbles. However, when bubbles collide, their spherical symmetry is broken and gravitational waves are generated. The corresponding spectrum can be computed using the envelope approximation [103]. In this approximation, the GWs are sourced from the fraction κ φ of the latent heat stored in thin shells around the bubble walls, taking into account only the uncollided part of the shells, i.e. the envelope of the collided bubbles.
Based on numerical simulation, the corresponding GW spectrum is given by [100]  and κ φ = ρ φ /ρ vac is the efficiency factor for the conversion of the vacuum energy density into scalar gradient energy. The scalar contribution can typically be neglected, unless the transition proceeds in the run-away regime. Here, due to the friction exerted by the lepton number gauge boson the bubbles do not run away [99]. We thus set κ φ = 0 and neglect the scalar contribution in the following. The second source comes from acoustic production. Since the scalar bubbles expand in the primordial plasma of the early Universe, their expansion excites the particles in the plasma and induces sound waves. If the wall velocity is less than the speed of sound in the fluid, c s = 1/ √ 3, the sound shell forms a deflagration in front of the wall, whereas for supersonic bubble walls, the shell becomes a detonation behind the wall. Similar to the scalar bubbles, the sound shells radiate gravitational waves upon collision. However, in contrast to the former case, it is not sufficient to consider only the overlap of collided shells. The collided parts of the shells do not disappear, but pass through one another, giving rise to a long-lasting (compared to the scalar bubble collisions) source of GW with a lifetime on the order of a Hubble time [101]. As a consequence, the sound-wave contribution typically dominates the GW spectrum.
Again, the spectrum is obtained by fits to numerical simulations, giving [101] with the spectral shape The efficiency factors for the conversion of the vacuum energy density into bulk motion for relativistic bubbles in the non-runaway case is given by [104] κ v α 0.73 + 0.083 √ α + α −1 . (6.10) Last but not least, GWs are also generated by turbulences. The expanding bubbles can induce eddies in the fluid. Since the plasma is ionized, these give rise to magnetohydrodynamic (MHD) turbulences. However, simulations indicate that turbulence effects are negligible compared to the sound waves [101].
Assuming Kolmogorov-type turbulences, the turbulently generated spectrum has the form [102,105]  where with the red-shifted Hubble rate at GW generation h * = 16.5 µHz T * 100 GeV g * 100 For the efficiency factor for the conversion into MHD turbulence, simulations suggest κ turb = κ v with ∼ 5 − 10%. We here follow [88] and take = 0.05. An example of the GW spectrum generated by the lepton number breaking phase transition for a scalar mass of m φ = 200 GeV and a gauge boson mass of m Z = 1.4 TeV (with v Φ = 2 TeV and L = −1/2, neglecting the heavy lepton contributions) is shown in figure 19 (blue curve). The contributions from acoustic production (green) and MHD turbulence (red) are indicated by dashed lines. For this choice of parameters, the transition occurs at a nucleation temperature of T n ∼ 200 GeV with a peak frequency of the spectrum of f ∼ 22 mHz. The spectra for further benchmark points and the corresponding values of T c , T n , α, and β/H can be found in appendix C.

Detectability
The stochastic GW background generated by a first-order PT can in principle be detected by GW observatories. Whereas current earth-based experiments such as LIGO and Virgo

JHEP02(2019)048
are not sufficiently sensitive to constrain first-order PTs in the frequency range we are interested in (due to seismic noise), space-based interferometers provide promising prospects for the detection of a stochastic background. We therefore evaluate detection prospects for LISA [94], B-DECIGO [106], DECIGO [107] and BBO [108].
The Laser Interferometer Space Antenna (LISA) [94] is an upcoming space-based GW interferometer consisting of six laser links on three satellites arranged in a regular triangle with an arm length of 2.5 × 10 9 m and a mission life-time of 4 years. Its launch is scheduled for 2034 [97]. The Big Bang Observer (BBO) [108] is a proposed follow-up experiment consisting of four LISA-like detectors, again arranged in a regular triangle. The fourth detector is co-located with one of the other detectors, forming a star-of-David configuration. A similar configuration has been proposed for the Japanese Deci-Hertz Interferometer Gravitational Wave Observatory (DECIGO) [107], which is planned to be launched in the 2030's [109]. Its scaled-down version B-DECIGO (Basic/Base-DECIGO) [106] supposedly starts in the late 2020's [109].
The detectability of a stochastic background by a given experiment can be evaluated by calculating the signal-to-noise ratio (SNR) ρ. If the SNR is greater than a threshold value ρ thr , we consider the background as detectable.
The SNR can be calculated from the GW spectrum and the experiment's noise spectrum. For a single detector such as LISA and B-DECIGO, it is given by [88,110] whereas for a network of detectors like BBO or DECIGO it reads [97,110] Here, T is the mission life-time in seconds, (f min , f max ) is frequency range of the experiment, and h 2 Ω n is the experiment's (effective) noise. For the noise curves we use the expressions given in [111] (LISA), [112] (BBO), [113] (DECIGO), and [106] (B-DECIGO). We further assume T = 4 yrs and ρ thr = 10 [88] (ρ thr = 8 [106] for B-DECIGO).
The sensitivity curves of the experiments are shown as gray lines in figure 19. Note that these curves are not the noise curves, but the power-law integrated curves [110]. These indicate that a GW background should be detectable by the experiment if the spectrum touches or reaches into the region above the respective sensitivity curve. Consequently, the spectrum shown in the figure can be probed by all four experiments. Figure 20 shows the parameters α and β/H * for the lepton number breaking phase transition as a function of the φ and Z mass in the simplified case of negligible portal coupling λ p considered in section 5.2, calculated using CosmoTransitions [93]. Most choices of masses give rise to a rather short first-order phase transition (high β/H * ) with few energy released (low α). However, large values of α and small values of β/H * can be obtained in the m Z 2m φ region, which is the region we identified to give rise to strong first-order PTs in section 5.2. As the amplitude of the sound-wave contribution to the GW spectrum (6.7) is proportional to α 2 /(1 + α) 2 and H * /β, this is indeed the region in which the stochastic background can be expected to be detectable.

Results
The corresponding parameter points for which the stochastic GW background generated by the lepton number breaking PT is accessible to space-based GW interferometers are depicted in figure 21(a). The red, purple, blue, and green regions can be detected by LISA, B-DECIGO, DECIGO, and BBO, respectively, whereas in the gray region the generated GW background in not detectable. If a parameter point is detectable by more than one experiment, the color corresponds to the experiment named first in the list above.
If the φ Yukawa couplings (and the Higgs portal coupling) are neglected, LISA and B-DECIGO are only sensitive in the m Z m φ margin of the parameter space that has a first-order PT. DECIGO can probe a small portion of the parameter space, also only close to the m Z m φ edge, BBO is slightly more sensitive. Still, the majority of the parameter space is inaccessible to GW experiments. vanishing Yukawa couplings. B-DECIGO can probe additional parameter points, mostly for low scalar masses. DECIGO and BBO can reach m Z 1 TeV for m φ ∼ 50 GeV and m Z 2.8 TeV for m φ ∼ 1 TeV. Last but not least, figure 21(d) shows the detectability of the GW from the lepton number breaking PT requiring that the DM accounts for the full thermal abundance measured by Planck, assuming m HL = 1.5 × m DM as in figure 18. Again, the effects of the dark leptons significantly enhance the parameter space to which future space-based GW observatories are sensitive. The dashed lines indicate the exclusion reach of XENON1T (cf. figure 5(b)) for DM mixing angles of sin θ DM = 0.015, 0.02, 0.22 and 0.024. The white region in the upper left part of the plot is excluded as the lepton number gauge group remains unbroken. Although not specifically mentioned, this applies to all sub-figures of figure 21.

Summary
In this work we have studied an extension of the SM in which lepton number is promoted to a U(1) gauge group. We updated the collider and dark matter limits in [4] and added further constraints. In particular, we considered the lepton number breaking phase transition. We investigated the parameter region in which the transition is of first order and evaluated the detectability of the generated stochastic background of gravitational waves.
Beside the lepton number gauge boson and the scalar field that spontaneously breaks U(1) , the model features additional leptonic states. Their presence is forced upon us by the necessity to cancel the gauge anomalies associated with U(1) . The lightest of these additional dark leptons is stable, unless specific values of the lepton number charges are assumed. Hence, the model naturally provides a dark matter candidate.
Assuming that the DM is a thermal relic, we identified the regions of the parameter space in which the DM candidate can account for the full abundance measured by Planck. We found that the correct relic density can be reproduced for a broad extent of DM masses in the O(100 GeV) to TeV range. This typically requires choosing m Z ∼ 2m DM .
Direct and indirect detection experiments put limits on the DM interactions with SM fields. Direct detection constrains the various mixings that can give rise to DM-quark interactions. These are the SM doublet admixture into the singlet DM characterized by θ DM , the kinetic mixing parameter of the lepton number and hypercharge gauge groups, and the mixing angle θ H between the SM and the lepton number Higgs. XENON1T can exclude and θ DM in the percent range, and sin θ H ∼ O(0.1); LZ and DARWIN can improve the limits by roughly an order of magnitude. Indirect detection on the other hand mainly probes Z mediated DM-SM interactions. However, even the next-generation CTA is only sensitive for lepton number charges as large as L = 3/2.
We also investigated collider limits. The most important ones are LEP limits. These put a lower bound on the lepton number breaking scalar VEV v Φ 1. Finally, we investigated the lepton number breaking phase transition in the early Universe. If the portal coupling between the SM and dark Higgs is sufficiently small, the lepton number and electroweak phase transitions happen independently from one another. Due to the VEV hierarchy imposed by the LEP constraints, the former typically occurs first. We found that in a large fraction of the parameter space the lepton number transition is first order. It can thus generate a stochastic background of gravitational waves. We calculated the corresponding GW spectrum and evaluated the detection prospects for future space-based GW observatories.
Whereas LISA can only probe a rather small fraction of the parameter space, its possible successors BBO and DECIGO are able to explore a significant fraction of the parameter points that give rise to a first-order PT. Notably, the exotic leptons significantly enhance the detection prospects, particularly when requiring that the measured relic abundance is reproduced. This is due to two effects. First, the presence of additional particles lowers the nucleation temperture, and second, the fermionic contributions restore the broken minimum in a part of the parameter space in which the bosonic corrections alone would shift the vacuum back to the origin.
Further interesting effects may arise if one considers non-vanishing portal couplings between the dark and SM scalar sectors. The transition can then proceed diagonally in field space, breaking the electroweak and lepton number gauge groups simultaneously. We leave this subject for future work. Another possible direction would be the investigation of the phase transition in the context of Baryogenesis. Also, a more elaborate evaluation of exclusion or discovery prospects at future colliders such as the ILC or a 100 TeV pp collider is needed.

A Goldstone divergences
In this appendix we address the cancellation of the IR divergence in the second derivative of the effective potential, originating from the vanishing Goldstone mass in Landau gauge. We follow the treatment in [92]. This leads to the renormalization condition (5.15). We calculate the self-energy Σ(p 2 ) of the lepton number breaking scalar φ in Landau gauge, using dimensional regularization in D = 4−2 dimensions. More precisely, we are interested in the difference of the self-energy evaluated at the scalar mass p 2 = m 2 φ and at p 2 = 0, ∆Σ ≡ Σ(m 2 φ ) − Σ(0), where p 2 is the external momentum squared. Figure 22. Tadpole diagrams contributing to the self-energy of φ.

A.1 Scalar self-energy
We consider the Lagrangian We rewrite the complex scalar Φ in terms of its real and imaginary parts and its vacuum expectation value v Φ as Φ = 1 √ 2 v Φ + φ + iω 0 . The (bare) self-energy of φ receives contributions from loops of φ itself, the Goldstone boson ω 0 , and the gauge boson Z . As the tadpole diagrams depicted in figure 22 are independent of the external momentum, we only need to evaluate the remaining bubble diagrams. These are We perform the Passarino-Veltman reduction [116] of the corresponding integrals using FeynCalc 9.3.0 [117,118], yielding where λ S = 6λ Φ for S = φ and λ S = 2λ Φ for S = ω 0 , respectively. The scalar one-and two-point functions A 0 and B 0 are defined by The renormalized self-energy is related to the bare one by Σ R p 2 = Σ 0 p 2 + δm 2 − p 2 δZ , (A.9) where δm 2 and δZ are the mass and field renormalization counter-terms for φ, L ⊃ δZ∂ µ φ∂ µ φ − δm 2 φ 2 . (A.10) When calculating ∆Σ, δm 2 cancels in the difference but δZ remains, i.e.
In particular, this means that ∆Σ is independent of the renormalization conditions we impose on the counter-terms δm 2 and δλ when calculating the effective potential. We can now fix δZ by requiring canonical normalization of the field φ, i.e.
Again, the tadpole diagrams in figure 22 do not contribute as they are independent of p 2 . Finally, the difference in the self-energy used in the renormalization condition (5.15) is obtained by plugging (A.8) and (A.12) into (A.11). We use LoopTools 2.13 [119,120] to evaluate the finite part of the scalar integrals and their derivatives.

A.2 On-shell renormalization of the effective potential
The momentum-dependent mass of φ is given by where m φ,R is the renormalized mass parameter in the Lagrangian, which is related to the physical (pole) mass m 2 φ ≡ m 2 φ (m 2 φ ) by (A.14) Since the effective potential is defined at vanishing external momentum, we now impose the conditions We now further want the VEV v Φ and the scalar mass m φ to be identical to the values inferred from the tree-level potential, i.e. contribution (equations (6.7) and (6.11)) the amplitudes of the spectra are proportional to v w and the peak frequencies shift as 1/v w , i.e. order one changes in the wall velocity only have an order one effect on the spectrum and peak frequencies. Hence, the detectability of the generated stochastic background eventually only has a mild dependence on v w . As a detectable signal further typically requires strong transitions with large wall velocities, we simply take the most optimistic estimate v w = 1. Figure 23 shows the dependence of the detectability on the bubble wall velocity in the case of negligible portal and Yukawa couplings. Compared to figure 21(a) above where v w = 1 was assumed, here we show the sensitivity for a slightly lower wall velocity (v w = 0.9), a wall velocity close to the speed of sound (v w = 0.6), and subsonic bubbles (v w = 0.3 and v w = 0.1). Again, the colored regions are detectable by LISA (red),  Table 3. Benchmark values of the scalar and Z masses, along with the DM and exotic lepton masses for the benchmark points D1 -D6.
B-DECIGO (purple), DECIGO (blue) and BBO (green). The GW spectra generated by the first-order phase transitions in the gray region are not detectable. For supersonic bubbles, the detectable regions barely change when varying v w . Taking v w to subsonic values decreases the sensitivity visibly.

C GW benchmark points
We here present results for 24 benchmark points. We consider the 6 combinations of scalar and Z masses listed in table 3. At each mass combination we take 4 different scenarios for the exotic leptons. Scenario A neglects the fermionic contributions, whereas B and C assume m DM = 200 GeV, m HL = 210 GeV and m DM = 500 GeV, m HL = 1 TeV, respectively. Again, m HL ≡ m e 4 = m e 5 . Finally, scenario C takes the values listed in table 3, for which the DM accounts for the full relic abundance measured by Planck. All benchmarks assume v Φ = 2 TeV and L = −1/2. All mixing angles are set to zero. We have computed the critical and nucleation temperatures T c and T n as well as the transition parameters α and β for the considered benchmark points. Table 4 lists the corresponding results obtained using CosmoTransitions [93]. The respective GW spectra are depicted in figure 24, assuming a bubble wall velocity of v w = 1. The sensitivities of LISA, B-DECIGO, DECIGO, and BBO are also shown, with the SNRs included in  Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.