Gravitational Wave Signals of Pseudo-Goldstone Dark Matter in the $\mathbb{Z}_{3}$ Complex Singlet Model

We study pseudo-Goldstone dark matter in the $\mathbb{Z}_{3}$ complex scalar singlet model. Because the direct detection spin-independent cross section is suppressed, such dark matter is allowed in a large mass range. Unlike in the original model stabilized by a parity, due to the cubic coupling of the singlet the $\mathbb{Z}_{3}$ model can accommodate first-order phase transitions that give rise to a stochastic gravitational wave signal potentially observable in future space-based detectors.


I. INTRODUCTION
One of the best candidates for the dark matter (DM) is a scalar singlet [1,2]. The properties of singlet DM have been studied in detail [3][4][5][6] (see [7,8] for recent reviews; see also Refs. therein). The non-observation of dark matter by direct detection experiments [9][10][11], however, puts severe bounds on models of DM comprised of weakly interacting massive particles (WIMP), pushing the mass of singlet scalar DM -except around the Higgs resonance -over 1 TeV.
A way to suppress the direct detection cross section is to consider as the DM candidate a pseudo-Goldstone with derivative couplings to CP-even states [12] (the DM phenomenology of the imaginary part of the complex scalar singlet was first considered in [3,13], but only in the Higgs resonance region). Because the velocity of DM particles in the Galaxy is small, the prospective signal is suppressed by vanishing momentum transfer. This result remains practically unaffected when loop corrections to the direct detection cross section are taken into account [14,15]. Despite that, it is possible for pseudo-Goldstone DM to show up at the LHC [16][17][18][19].
In this class of models, the global U (1) symmetry is explicitly, but only softly broken into a discrete subgroup which is then broken spontaneously. In Ref. [12], the U (1) group was explicitly broken into Z 2 symmetry. We study the consequences of breaking U (1) into Z 3 . The model admits two phases that produce a dark matter candidate. Firstly, the unbroken Z 3 symmetry stabilizes S as a dark matter candidate. Secondly, with the broken Z 3 symmetry, the imaginary part of S, denoted by χ, is still stable due to the S → S † symmetry of the Lagrangian. The main difference between the Z 2 and Z 3 pseudo-Goldstone DM models is that the potential of the latter contains a cubic S 3 term. There are other possible cubic terms [20][21][22], but they do not respect the Z 3 symmetry.
To this date, the Z 3 -symmetric complex singlet model has been studied in the unbroken phase. Originally, the * kristjan.kannike@cern.ch † kaius.loos@gmail.com ‡ martti.raidal@cern.ch model was proposed in the context of neutrino physics [23]. Detailed analysis of DM phenomenology was carried out in Ref. [24]. Indirect detection of Z 3 DM was considered in [25,26]. The cubic coupling can contribute to 3 → 2 scattering for Z 3 strongly interacting (SIMP) DM [27][28][29][30]. The effects of early kinetic decoupling were studied in [31]. The Z 3 symmetry has been considered as the remnant of a dark U (1) local [32][33][34] or global [35] symmetry. In the unbroken phase, the direct detection cross section can also be suppressed if the DM relic density is determined by semi-annihilation processes [36][37][38][39][40][41]. However, the suppression is not as large as for pseudo-Goldstone DM. The discovery of gravitational waves (GWs) by the LIGO experiment [42,43] opened a new avenue to probe new physics. First-order phase transitions generate a stochastic GW background [44][45][46] which may be discoverable in future space-based GW interferometers [47,48]. While the SM Higgs phase transition is of second order [49,50] and does not generate a GW signal, in models with extended scalar sector the first-order phase transition in the early Universe can become testable by observations. For a recent review on phase transitions and GWs, see Ref. [51].
GWs from beyond-the-SM physics with a scalar singlet have been studied in detail. These models admit a two-step phase transition that can be of the first order [52][53][54][55][56][57][58][59] and can potentially produce a measurable GW signal [60][61][62][63][64][65][66][67][68][69][70][71]. However, in the Z 2 pseudo-Goldstone model, all phase transitions leading to the correct vacuum are of the second order [72], yielding no stochastic GW signal and excluding the additional potentially powerful experimental test of this class of the SM models. The phenomenology of phase transitions and GWs of the Z 3 complex singlet model was studied in detail in Ref. [73]. This work, however, disregarded the singlet as a DM candidate, because in the unbroken phase a sizeable GW signal is incompatible with correct relic density.
The goal of this paper is to study the nature of phase transitions and GW signals in the complex scalar singlet model in which a U (1) symmetry softly broken into its Z 3 subgroup. Just like in the Z 2 case, the elastic scattering cross section pseudo-Goldstone DM with matter can be small enough for the DM to be well-hidden from direct detection. We find, however, that unlike in the Z 2 case, strong first-order phase transitions can take place because the potential contains a cubic term. The resulting stochastic GW background is potentially discoverable in future space-based detectors such as LISA and BBO.
The paper is organized as follows. We introduce the model in Sec. II. Various theoretical and experimental constraints are discussed in Sec. III. DM relic density, phase transitions and the predictions for the direct detection and GW signals are treated in Sec. IV. We conclude in Sec. V. The details on the effective potential and thermal corrections are relegated to Appendix A.

II. Z3 COMPLEX SINGLET MODEL
The most general renormalizable scalar potential of the Higgs doublet H and the complex singlet S, invariant under the Z 3 transformation H → H, S → e i2π/3 S, is given by where only the cubic µ 3 term softly breaks a global U (1) symmetry. Notice that the Z 3 symmetry precludes any quartic couplings that would result in a hard breaking of the U (1). The potential Eq. (1) -just like in the original Z 2 pseudo-Goldstone DM model -has an additional discrete Z 2 symmetry, S → S † .
In the unitary gauge, we parameterize the fields as The S → S † is then equivalent to χ → −χ, which makes χ stable even as the Z 3 symmetry is broken. The model thus admits two different DM candidates: in the unbroken Z 3 phase the complex singlet S is a DM candidate, while in the broken Z 3 phase it is its imaginary part, the pseudo-Goldstone χ, which can be the DM. We concentrate on the latter case. Without loss of generality, the vacuum expectation value (VEV) of S can be taken to be real and positive, because the degenerate vacua where χ has a non-zero VEV are related to the real vacuum by Z 3 transitions. This choice corresponds to a negative value of the parameter µ 3 .
The stationary point conditions are We see that the model admits four types of extrema: the origin as O ≡ (0, 0) is fully symmetric; the vacuum H ≡ (v h , 0), with Higgs VEV only, spontaneously breaks the electroweak symmetry; the vacuum S ≡ (0, v s ), with S VEV only, spontaneously breaks Z 3 ; our vacuum HS ≡ (v h , v s ) breaks both symmetries.
The mass matrix of CP-even scalars in our HS vacuum is given by The mass matrix (5) is diagonalised by an orthogonal matrix The mixing angle θ is given by The mixing of the CP-even states h and s will yield two CP-even mass eigenstates h 1 and h 2 The mass of the pseudoscalar χ is, taking into account the extremum conditions, which is proportional to µ 3 as it explicitly breaks the U (1) symmetry. We express the potential parameters in terms of physical quantities in the zero-temperature vacuum, that is the masses m 2 1 and m 2 1 of real scalars, their mixing angle θ, pseudoscalar mass m 2 χ , and the VEVs v h and v s : Both the Higgs doublet and the singlet will get a VEV, with the Higgs VEV given by v h = v = 246.22 GeV. We identify h 1 with the SM Higgs boson with mass m 1 = 125.09 GeV [74].

III. THEORETICAL AND EXPERIMENTAL CONSTRAINTS
We impose various theoretical and experimental constraints on the parameter space of the model. First of all, the potential (1) is bounded from below if Secondly, we require the couplings to be unitary and perturbative. The unitarity constraints in the s → ∞ limit are given by where the last condition, in the λ SH = 0 limit, yields |λ H | 4 3 π and |λ S | 2π. We also calculate unitarity constraints at finite energy with the help of the latest version [75] of the SARAH package [76][77][78][79]. Scattering at finite energy allows to set a bound on the cubic coupling µ 3 .
To ensure validity of perturbation theory, loop corrections to couplings should be smaller than their tree-level values. The model is perturbative [80] if |λ H | 2 3 π, |λ S | π and |λ SH | 4π.
We require that our HS vacuum be the global one: this implies that which for sin θ ≈ 0 is approximated by m χ 3m 2 . This bound appears to be stronger than the constraint on the cubic coupling from unitarity at finite-energy scattering.
If the mass of x ≡ χ or h 2 is less than m h /2, then the Higgs invisible decay width into this particle is given by with The invisible Higgs branching ratio is then given by which is constrained to be below about 0.24 at 95% confidence level [81,82] by direct measurements and below about 0.17 by statistical fits of all Higgs couplings [83,84]. In extended Higgs models [85] the mixing phenomenology can be more complicated.
The mixing angle between h and s is constrained from the measurements of the Higgs couplings at the LHC to | sin θ| ≤ 0.37 for m 2 m h and | sin θ| ≤ 0.5 for m 2 m h [86].
Last, but not least, we require that the relic density of χ be equal to the value Ωh 2 = 0.120 ± 0.001 from recent Planck data [87].

IV. DIRECT DETECTION AND GRAVITATIONAL WAVE SIGNALS
For the scan of the parameter space we choose m χ , m 2 and sin θ as the free parameters, while v s is used to fit the DM relic density. We generate the free parameter values in the following ranges: m χ ∈ [25, 1000] GeV, m 2 ∈ [25, 4000] GeV and sin θ ∈ [−0.5, 0.5] GeV. We use the micrOMEGAs package [88] to fit the DM relic density and CosmoTransitions [89] for the calculation of phase transitions and the Euclidean action to determine tunnelling rates. The stochastic gravitational wave signal was calculated using the formulae of Ref. [90] for the non-runaway case.
The micrOMEGAs is also used to compute predictions for direct detection signals, to which we add the tiny loop correction. Unlike in the Z 2 pseudo-Goldstone DM model, the tree-level direct detection DM amplitude contains a term that does not vanish at zero momentum transfer, so For a large part of the parameter space, however, this term is small enough, which allows one to explain the negative experimental results from DM direct detection experiments for a wide range of pseudo-Goldstone DM mass.
In the high temperature approximation (A8), the mass terms acquire thermal corrections: with the coefficients For numerical calculations of gravitational wave signals, exact expressions are used. We present in Fig. 1 the results of our scans for the direct detection signal. In the left panel we present for illustration the behavior of spin-independent direct detection cross section σ SI as a function of DM mass m χ for fixed values of m 2 = 300 GeV and sin θ = 0.05. Red lines show the current limit from the XENON1T experiment [10] and the future bound from the projected XENONnT experiment [91]. While the usual Z 2 scalar singlet DM is excluded below O(TeV) masses, except for the Higgs resonance region, the suppression of the pseudo-Goldstone cross section in this model allows for lighter DM candidates. The two pronounced dips in the cross section are given by the Higgs resonance at m h /2 and the h 2 resonance at m 2 /2. The direct detection signal varies over a wide range and is roughly proportional to the cubic coupling µ 3 . Because the amplitude is proportional to sin θ cos θ, larger mixing angles also produce larger direct detection signal. The right panel in Fig. 1 presents the results of the whole scan, coded in grey-scale by the mixing angle | sin θ|; empty points are excluded by the constraint on the Higgs invisible branching ratio. Note that the BR inv can exclude points with m χ > m h /2, because the Higgs boson can decay invisibly also to two h 2 if it is light enough. Early kinetic decoupling [92,93] may additionally enhance BR inv several times [31], but in practice this would not change the parameter space of GW signals. In particular, all our points with a potentially measurable GW signal have m χ , m 2 > m h /2. The upper bound on σ SI arises from the requirement the HS vacuum be global, which bound is stronger than that from unitarity that also constrains the points from above.
In Fig. 2 we present the predicted stochastic GW signal together with the predicted sensitivity curves of the future LISA and the BBO satellite experiments. To avoid cluttering the plot, we only show the peak power of each GW spectrum. In the left panel, light red points are excluded by the XENON1T direct detection constraints, while the darker green points are still allowed by direct detection; all points not satisfying the other constraints have been excluded. In the right panel, we depict the dependence of the GW signal on the cubic cou-pling µ 3 . The color code shows the nature of the phase transition that yields the strongest signal in a sequence of phase transitions. The largest signal is produced by the O → S transitions, enhanced by a sizable cubic coupling. Because m χ ∝ µ 3 , there is a similar dependence on m χ /m 2 . The strongest signals can be produced at roughly m χ ≈ 2.3m 2 .

V. CONCLUSIONS
We have studied the Z 3 complex scalar singlet DM model, where only the cubic coupling of the singlet explicitly breaks a global U (1) symmetry. The model has two phases with stable DM. In the phase where the Z 3 is spontaneously broken, the residual CP-like Z 2 symmetry stabilizes the imaginary part of the complex singlet S as a pseudo-Goldstone DM candidate. The DM direct detection cross section can be considerably suppressed by the small momentum transfer or the resonance at half the mass of the heavy singlet-like particle.
While we may be unable to discover the pseudo-Goldstone DM via direct detection, it may be testable by other means. In this model a stochastic gravitational wave background can arise from the first-order phase transitions due to the presence of a cubic coupling for the singlet. The strongest signals, which can potentially be observed by the future BBO experiment, are produced in the O → S phase transition (if present) at m χ ≈ 2. is smaller (at fixed cubic coupling).

ACKNOWLEDGEMENTS
We would like to thank Matti Heikinheimo and Christian Gross for useful comments. This work was supported by the Estonian Research Council grant PRG434, the grant IUT23-6 of the Estonian Ministry of Education and Research, and by the European Union through the ERDF Centre of Excellence program project TK133.

Appendix A: Thermal Effective Potential
At one-loop level, the quantum corrections to the scalar potential in the MS renormalisation scheme are given by where n i are the degrees of freedom of the i-th field, m i are field-dependent masses and the constants c i = 3 2 for scalars and fermions and c i = 5 6 for vector bosons. The masses and degrees of freedom n i of the fields are given in Table I. The field-dependent masses of h and s are given by the eigenvalues m 2 1,2 of the mass matrix of the CP-even eigenstates with elements given by We neglect the contributions of the Goldstone bosons G 0 and G ± . To calculate the effective potential in case of negative field-dependent masses, we substitute ln m 2 i → ln |m 2 i |, which is equivalent to analytical continuation [94]. We set the renormalisation scale to µ = M t .
We add a counter-term potential, δV = δµ 2 H |H| 2 + δλ H |H| 4 + δµ 2 S |S| 2 + δλ S |S| 4 + δλ SH |S| 2 |H| 2 + δµ 3 2 (S 3 + S †3 ) + δV 0 , in order to fix the VEVs and the mass matrix in our HS minimum to their tree-level values. The thermal corrections to the potential are given by where J ∓ = ± ∞ 0 dyy 2 ln 1 ∓ exp − y 2 + x 2 , (A7) with the − sign applied to bosons and the + sign to fermions. In the high-temperature limit m/T 1, the The full thermally corrected effective potential is then The counter-term potential δV is chosen such as to keep quantum corrections to the masses, to mixing between h and s and to the VEVs zero. The counter-terms are given by where we take h = v, s = v S , χ = 0 after taking the derivatives.