Pseudo-Goldstone dark matter: gravitational waves and direct-detection blind spots

Pseudo-Goldstone dark matter is a thermal relic with momentum-suppressed direct-detection cross section. We study the most general model of pseudo-Goldstone dark matter arising from the complex-singlet extension of the Standard Model. The new U(1) symmetry of the model is explicitly broken down to a CP-like symmetry stabilising dark matter. We study the interplay of direct-detection constraints with the strength of cosmic phase transitions and possible gravitational-wave signals. While large U(1)-breaking interactions can generate a large direct-detection cross section, there are blind spots where the cross section is suppressed. We find that sizeable cubic couplings can give rise to a first-order phase transition in the early universe. We show that there exist regions of the parameter space where the resulting gravitational-wave signal can be detected in future by the proposed Big Bang Observer detector.

The paper is organised as follows: We introduce the general complex-singlet model in section 2. In section 3 we discuss DM phenomenology including direct and indirect detection. Cosmic phase transitions are considered in section 4, and we conclude in section 5. The appendices contain some more technical details: In appendix A we give the formulae for shifting away the linear term in the potential, appendix B lists the annihilation cross sections, appendix C shows how the one-loop contribution to the direct-detection cross section is modified in the presence of cubic symmetry-breaking terms, and in appendix D we list the one-loop renormalisation group equations for the model.

General complex-singlet model
We consider a model where the scalar sector consists of the SM Higgs doublet, H, and a complex singlet, S. The model is by construction CP-conserving [30], i.e. invariant under the CP-like transformation S → S * . We write the potential as where is invariant under a global U(1) transformation S → e iφ S, while the remaining part explicitly breaks the U(1) symmetry: (2. 3) The minimal Z 2 -symmetric pseudo-Goldstone DM model contains only the U(1)-symmetric potential and the explicit symmetry-breaking µ 2 S mass term; the potential (2.1) is the most general setup. In the unitary gauge the scalar multiplets are parametrized as , and S = 1 √ 2 (w + s + iχ), (2.4) where v = 246. 22 GeV is the usual electroweak vacuum expectation value (vev), and w is the vev of the singlet scalar. The λ S term produces an independent contribution only to the χ 4 vertex; elsewhere it can be absorbed by a redefinition of the other couplings. Consequently, we will set λ S = 0 in what follows. Minimising the potential, the mass of the CP-odd field, χ, is calculated to be The real part of the singlet, s, mixes with the neutral component of the doublet, h, resulting in two mass eigenstates, h 1 and h 2 . The particle h 1 = h cos θ − s sin θ is identified as the SM-like Higgs boson with mass m 1 = 125.1 GeV, and the orthogonal linear combination h 2 is another CP-even scalar with mass m 2 . The mixing angle θ is given by tan 2θ = − 8(λ HS + λ HS )vw 2 + 4µ HS vw 8λ H v 2 w − 8(λ S + λ S )w 3 + 3 √ 2(µ 3 + µ 3 )w 2 + µ HS v 2 + 4µ 3 1 . (2.6) We replace the parameters µ 2 H , µ 2 S , µ 2 S , λ H , λ S and λ HS appearing in the potential with physical parameters m 1 , m 2 , m χ , θ, v and w. Fixing the value of the electroweak vev and the mass of the Higgs boson to the known values reduces the number of independent parameters by two. The independent parameters are then further constrained by collider searches and cosmological observations. Let us briefly discuss the phenomenological constraints which are relevant for our model.
Stability of the potential and unitarity. To guarantee a stable vacuum, the potential has to be bounded from below. This is in particular relevant for large field values, where we can neglect dimensionful terms. Imposing the co-positivity condition [31,32] on the matrix of quartic couplings requires where λ HR ≡ λ HS + λ HS + 2 λ H (λ S + λ S ) > 0, λ HI ≡ λ HS − λ HS + 2 λ H (λ S + λ S ) > 0.
In addition to this, we ensure that the point is the global minimum of the potential. Unitarity of the S matrix for two-to-two elastic scattering constrains the values of combinations of λ-parameters in the potential at asymptotically large center-of-mass energies [33,34]. In our numerical analysis with the SARAH package [35], we determine the eigenvalues Λ i of the scattering matrix, require |Λ i | ≤ 1/2 and implement the resulting constraints on the quartic couplings. In our model this condition translates into The three remaining eigenvalues Λ 1,2,3 of the scattering matrix are the solutions to a cubic equation and lie in the interval [36,37] and constrains the parameter space in each decay channel. In particular, the signal strengths constrain the value of the mixing angle to satisfy cos 2 θ ≥ 0.9 [39]. The latest measurement of the width of an SM-like Higgs boson gives Γ h 1 tot = 3.2 +2.8 −2.2 MeV, with 95% CL limit on Γ tot ≤ 9.16 MeV [40]. In our model, the total width of the SM-like Higgs boson can be modified -if kinematically allowed -through new decay channels h 1 → χχ , (2.12) where the partial decay width for a new decay channel h 1 → XX is: Current experimental values provided by the ATLAS and CMS experiments [41,42] on the invisible branching ratio are where h 1 → inv. represents the SM-like Higgs decay to the DM candidate, χ. We will use the conservative limit Br(h 1 → inv.) < 0.23 in the following analysis.
Relic density measurements. The obsevations by the Planck satellite [43] show the abundance of DM to be Ω c h 2 = 0.120 ± 0.001, (2.15) where h is the dimensionless Hubble parameter, H = 100 h km/s/Mpc. In terms of the relative DM abundance, f rel = Ω χ h 2 /(Ω c h 2 ), we take f rel = 1 meaning that χ constitutes all of the expected DM relic density.
Taking into account these preliminary constraints, we then consider direct and indirect detection of the DM candidate in our model to identify regions of the parameter space surviving all the above constraints.

Direct and indirect detection
Let us now consider the cross sections relevant for the phenomenology of the model. Our focus will be on understanding the suppression of the direct-detection cross section, and how the situation changes in the presence of various symmetry-breaking terms.

Tree-level cross section
Let us first review the suppression of the direct-detection cross section at tree level. The CP-even scalar mass eigenstates couple to χ as Here the couplings are and we have defined and The two CP-even scalar mass eigenstates couple to the nucleon, N , via the Higgs-boson Yukawa couplings as GeV is the nucleon mass, and we use f N = 0.3 [44][45][46] for the effective Higgs-nucleon coupling. The spin-independent direct-detection cross section is given by where we have defined the effective DM-nucleon coupling λ SI as Writing the scalar couplings explicitly, cf. eq. (3.2), the effective direct-detection coupling in the t → 0 limit becomes As discussed in the introduction, the tree-level direct-detection cross section vanishes in the t → 0 limit for U(1)-invariant interactions or if only the Z 2 -symmetric µ 2 S mass term is included. As explicitly shown by the above equation, this is no longer true if any other symmetry-breaking interaction terms in the potential, eq. (2.3), are present.
To extract the effect of the symmetry-breaking terms on the direct-detection cross section, it is instructive to study the DM-nucleus interaction regardless of the relic-density contribution of the DM candidate, χ. In figure 1, we show the spin-independent DMnucleus interaction cross section for the minimal Z 2 -symmetric model enhanced with only one non-zero symmetry-breaking term in each plot. Figure 2 shows regions allowed by the XENON1T bound [1,2] regardless of the relic-density contribution of χ for typical values of w = 250 GeV, m 2 = 200 GeV and sin θ = 0.1. In the left panel, all terms with odd powers of the singlet S in the symmetry breaking potential in eq. (2.3) have been set to zero, i.e. µ HS = µ 3 = µ 3 = µ 1 = 0, and in the right panel all doublet-singlet mixing terms in eq. (2.3) have been set to zero, i.e. µ HS = 0 and λ HS = 0.

Cancellation regions
Let us study the general formula for the effective direct-detection coupling in the t → 0 limit in eq. (3.8) in more detail. As noted before, this result is generally non-zero in the presence of any of the symmetry-breaking interactions in eq. (2.3). However, there are specific combinations of the symmetry-breaking parameters which lead to a suppressed λ SI in the t → 0 limit, thereby mimicking the behaviour of the minimal Z 2 -symmetric model. We shall now explore these cancellation conditions more closely.
Recall first, that in the non-linear representation, if only the U(1)-breaking µ 2 S mass term in eq. (2.3) is present, we have yielding the following effective coupling for direct detection, where p 1 and k 1 are the mo-  menta of the incoming and outgoing χ-particles, resp., which explicitly shows that the direct-detection cross section vanishes in the t → 0 limit. Let us now study how the cubics µ 3 and µ 3 affect the cross section and pseudo-Goldstone mass. For simplicity, we take here λ S = 0 and µ 3 1 = 0 and also set all symmetrybreaking terms involving the Higgs boson to zero. Representing the singlet field, S, as the cubic terms in the potential (2.1) can be written as Minimization of the potential yields  Figure 3: The potential, V 3 , and its first two derivatives with respect to χ at s = 0 for µ 3 = −9µ 3 . and the mass of χ is calculated to be (3.14) Let us now try to understand the origin of the contributions to the direct-detection cross section by explicitly relating the derivation of the Z 2 -symmetric case, eq. (3.10). The relevant interaction terms in our simplified case are given by where we used eq. (2.5) to obtain the last equality. The second term in the parenthesis of the last line of eq. (3.15) does not cancel by the mass coming from the derivative terms, and thus there is a contribution to direct detection from the cubic terms. However, note how the mass m χ , given by eq. (3.14) in the simplified case at hand, is proportional to the same combination of parameters as in the direct-detection cross section. This shows explicitly that also in this case the suppression of the direct-detection rate is tied to the pseudo-Goldstone nature of the χ field.
In the special case µ 3 = −9µ 3 , when the tree-level direct-detection cross section cancels, the contribution to the m χ from the cubic terms goes to zero as well, thereby implying a more symmetric vacuum than just an accidental cancellation. This can be understood by the form of the χ-potential at the vacuum, eq. (3.12). To illustrate this explicitly, we show the cubic part of the potential and its derivatives for the special parameter domain µ 3 = −9µ 3 in figure 3. As we noted above, the second derivative of V 3 with respect to χ vanishes when evaluated in the vacuum. However, notice that the cubic contributions along the s direction are not zero at χ = 0 even for µ 3 = −9µ 3 implying that the full Lagrangian does not have an enhanced symmetry in this limit.
In the general case of eq. (3.8), setting the combinations of parameters which appear in the direct-detection cross section, to zero leads to suppression of the direct-detection rate. It can be shown that these same combinations also appear in m χ . To illustrate these conclusions, we show in the left panel of figure 4 the contours of the ratio of tree-level σ SI to the XENON1T upper limit on the DM-nucleon cross section for the parameter combination 9µ 3 + µ 3 . Thus, contours with values of one or less are allowed. For this plot, we have chosen the typical values of w = 250 GeV, m 2 = 200 GeV and sin θ = 0.1, while all other symmetry-breaking terms (except for µ 2 S ) are set to zero.

Direct-detection cross section at one loop
As shown in [28], in the case of the simplest U(1)-invariant model broken by only the Z 2symmetric mass term, the non-vanishing corrections to the direct-detection cross section in the t → 0 limit arise at one-loop order. When more general U(1)-breaking interactions are considered, they will yield O(t 0 ) contributions to the direct-detection cross section already at tree level. The allowed magnitude of these interactions at tree level is expected to be roughly similar to the size of the loop corrections due to the mass term [6]. The loop corrections arising from the symmetry-breaking interactions are then expected to be negligible. However, as we have discussed above, the contributions from the symmetry-breaking interactions are suppressed at tree level in specific parts of the parameter space. In such case the effect of loop corrections becomes again relevant. Therefore, we will briefly discuss the one-loop contribution in the presence of these cubic terms extending the analysis of ref. [28].
We will concentrate here on the case where the only non-zero cubics are µ 3 and µ 3 . This choice is motivated by simplicity, but also because we want to, in particular, study if the suppression of the direct-detection cross section for the specific choice µ 3 = −9µ 3 is preserved at the loop level.
The spin-independent cross section with the one-loop corrections at t → 0 limit can be written as where now λ tree SI is given by eq. (3.8) with µ 1 = µ HS = 0 GeV, and λ HS = λ S = 0. In general, λ 1L SI has a complicated analytic expression involving several loop functions. In the particular case δ ≡ 9µ 3 + µ 3 = 0 the expression simplifies significantly, and we show the result for illustration in appendix C. In the numerical computation we keep the full expression including the deviation from the δ = 0 limit.
While at tree level, the value of only the combination δ is relevant, at loop level also the individual values of the coupligns µ 3 and µ 3 become important, since the limit δ = 0 does not correspond to a symmetry at the Lagrangian level. We illustrate this in figure 4, where we plot the spin-independent cross section (again regardless of the relic-density contribution) at tree and one-loop levels as a function of m χ and δ for two representative values of µ 3 = 50, 500 GeV. For the numerical evaluation of the various loop functions, we use the pySecDec toolbox [47,48] with FORM optimization [49][50][51] and CUBA library for multi-dimensional integration [52,53]. Figure 4 shows that the suppression of the direct-detection cross section persists also at one-loop level for moderate values of the couplings µ 3 , µ 3 , but at larger values the interference with the symmetry-breaking mass term increases the one-loop cross section significantly at large m χ . Note however, that near the δ = 0 region for δ < 0 and large m χ , the tree-level and one-loop contributions interfere destructively, allowing for the region of suppressed direct-detection cross section for larger DM masses, too.
We conclude that the situation is analogous to the simplest Z 2 -symmetric model of pseudo-Goldstone DM: Starting from a model in the cancellation region of U(1)-breaking parameters discussed in the previous section, we expect deviation from the symmetric t → 0 limit of the direct-detection cross section by contributions of the order of the loop corrections discussed in this section.
Another feature arising from loop corrections is that the parameter combination δ = 9µ 3 + µ 3 may not stay zero under running of couplings. The β-functions of the model are given in appendix D. As a simple example, let us consider the case where all other symmetry-breaking interactions are set to zero except µ 3 and µ 3 . Then we have which explicitly shows that the running of δ is not multiplicative. Generally, the model should be viewed as a low-energy effective theory with the coefficients of the symmetrybreaking operators taking non-zero values constrained to be compatible with experiments and observations. Nevertheless, the cancellation regions we have discussed here may be interesting towards more complete model building, in particular with the relatively large regions of the parameter space near the δ = 0 limit with suppressed direct-detection cross section persisting even at one-loop level.

Indirect detection
Finally, to relate the present analysis to the results of ref. [6], we discuss the implications of indirect detection of DM for our model framework. The relevant indirect-detection constraints arise due to annihilation of DM in the dwarf galaxies that orbit our Milky Way. In these structures the DM is cold and therefore DM annihilations take place essentially at vanishing momentum: s → 4m 2 χ . The annihilation cross section to all final states, ff , W + W − , ZZ, h 1 h 1 , h 1 h 2 and h 2 h 2 is non-vanishing in this limit (unless kinematically forbidden). The annihilation cross sections for these processes are given in appendix B.
The constraints from the Fermi-LAT dwarf galaxy observations, based on the recent analysis [54], are presented in figure 5. The red curve shows where the model reproduces the observed relic density, f rel = 1. The excluded regions, shown in blue, are produced by comparing the annihilation cross section v rel σ χχ→ bb to the reported bb exclusion limit for m χ < m W ± . In the region m χ > m W ± , the dominant annihilation channel is to W + W − and in this region we compare the total DM annihilation cross sectionthat is, the combined cross sections χχ → W + W − , ZZ, h 1 h 1 , h 1 h 2 , h 2 h 2 -to the W + W − bound. The kinks in the plot occur at kinematic thresholds for each channel, that is at m χ = m W ± , m Z , m 1 , (m 1 + m 2 )/2, m 2 . We have set all the explicit U(1)-breaking terms to zero in figure 5, except for the mass term µ 2 S and the cubics µ 3 and µ 3 satisfying δ = 9µ 3 + µ 3 = 0.
The gamma ray flux originating from DM annihilations in a dwarf galaxy depends on the density profile of the DM halo. This effect is described via the so-called J-factor.
In the present analysis we use the constraints on χχ → bb and χχ → W + W − cross sections based on J-factors obtained in [54]. These updated bounds are weaker than the ones used in our previous work [6], 1 and therefore we find the model less constrained by the indirect-detection data. We conclude that, taking into account the uncertainties in the determination of the J-factors, the model is not presently constrained by indirect detection in the m χ > m 1 /2 region.

Thermal potential
So far we have seen how the pseudo-Goldstone DM model can account for the observed relic abundance, and how it is constrained by direct-detection experiments. Since the model consists of an extended scalar sector, it is natural to explore the finite-temperature phase transitions in the early universe and their phenomenological consequences, in particular for the GW signals relevant for the LISA, BBO or DECIGO satellites.
It turns out that the phenomenologically viable scenario is a two-step transition starting from the high-temperature vacuum in the singlet direction (h, s) = (0, w 0 ). With the linear and the cubic terms present -in the absence of an additional Z 2 symmetry -the symmetry in the s direction is not necessarily restored, and w 0 can be non-zero. In the first step, as temperature is lowered, another minimum forms in the singlet direction, (h, s) = (0, w 1 ), and the transition (0, w 0 ) → (0, w 1 ) occurs at the critical temperature T c . This transition is potentially of first order, and can produce GW signals to be searched for in the future space-based missions. The electroweak transition from (0, w 1 ) to (v 2 , w 2 ) happens at a significantly lower temperature T EWPT c T c , and finally evolves to the zerotemperature global minimum (v, w). In the phenomenologically viable parameter space for DM, this second transition is predominantly of second order, and does not produce detectable GW signals. In the following, we will concentrate on the former first-order transition.
Let us now turn to the quantitative analysis of this scenario. The effective one-loop potential reads where the tree-level potential is given by eq. (2.1). The second term, V 0 CW , is the T = 0 Coleman-Weinberg potential in the MS scheme, where g i denotes the number of degrees of freedom, M i (h, s) are the field-dependent masses, and µ 0 is the renormalisation scale (which we fix to be µ 0 = v). In this expression, F = 1 for fermions and 0 for bosons, C i = 3/2 for scalars, fermions and longitudinal polarizations of gauge bosons and 1/2 for transverse polarizations of gauge bosons. The one-loop finite-temperature corrections are given by where the upper signs correspond to bosons and lower signs to fermions. The last term in eq. (4.1) contains the finite parts of the counter terms that are fixed such that the scalar vevs and masses remain at their tree-level values at the minimum, eq. (2.8): such that the following renormalization conditions are satisfied: Finally, we use the thermally improved finite-temperature potential, which is obtained by adding to the field-dependent masses in eqs. (4.2), (4.3) the leading thermal corrections (see ref. [58] for a recent discussion): where the coefficients c i are given by c h =(9g 2 + 3g 2 + 12y 2 t + 24λ H + 4λ HS )/48, c s =(2λ HS + 2λ HS + 4λ S )/12, (4.7) c χ =(2λ HS − 2λ HS + 4λ S )/12.

Gravitational-wave signal and peak-integrated sensitivity curves
During a first-order phase transition, stochastic GWs are produced via three independent mechanisms: collisions of bubbles (b), sound waves in the plasma (s), and turbulence in the plasma (t). The resulting GW spectra can be approximately written in terms of a peak amplitude, Ω peak i , and a spectral shape, S i which depends on the peak frequency, f i , The peak amplitudes and peak frequencies depend on the characteristics of the phase transition which can be quantified in terms of the nucleation temperature, T n , the amount of energy density released relative to the radiation energy density [59] characterizing the strength of the transition, and in terms of which gives approximately the inverse duration of the transition [60]. Here, S is the Euclidean action of the bubble solution. The latent heat released during the phase transition is converted with efficiency κ b into the kinetic energy of the expanding bubbles, with κ s into the sound waves, and with κ t into the turbulent motion in the plasma. There is still no consensus in the literature on how latent heat released into the kinetic energy of the plasma, 1 − κ b , is subsequently transformed into sound waves and turbulence [27,59,[61][62][63]; here we choose the commonly used estimate for the turbulence fraction κ t = 0.1, and we fix the bubble wall velocity v w = 0.9. The peak amplitudes and peak frequencies depend on parameters v w and κ i and on the values of functions α and β evaluated at the nucleation temperature; explicit formulas for these in the runaway-bubble-in-plasma scenario are given in refs. [59,62] along with the spectral shape functions, S i , in terms of the peak frequencies.
The resulting GW signal needs to be compared with the noise spectrum of the experiment under consideration to obtain the signal-to-noise ratio (SNR) [64,65] where n det = 1 if the experiment consists of just one detector allowing only auto-correlation measurement, while n det = 2 if it is possible to cross-correlate signals of a detector pair. In the following, we will consider three satellite-borne GW interferometers: LISA [16,17], DECIGO [18][19][20][21], and BBO [22][23][24] with the planned configurations n det = 1 for LISA and n det = 2 for DECIGO and BBO. The noise spectra of these experiments are discussed in detail in ref. [66].
For the representation of the model parameter points in the GW signal region and the experimental reach of the above interferometers, we adopt the approach of peak-integrated sensitivity (PIS) curves put forward recently in refs. [66,67]. The advantage of this approach with respect to the conventional power-law-integrated sensitivity curves [68] is that it allows to represent each parameter point as a single point in the GW signal region and is thus well-suited for the purpose of a general scan of the parameter space to be discussed in the following section.
The key observation is that if the shape of the expected signal is known, as is the case of first-order phase transitions, the integration over the spectral shape can be carried out leaving SNR uniquely determined by the peak energy densities and peak frequencies that depend on the model-specific phase-transition quantities and no longer on the GW frequency. More specifically, the SNR in eq. (4.11) can be rewritten as where the integration over the frequency range has already been carried out implicitly: where i, j ∈ {b, s, t} and the mixed peak amplitudes are defined as geometric means, (4.14) For details of the approach, see refs. [66,67].

Scan of the parameter space and results
We scan the parameter space with non-zero soft-breaking mass and cubic interactions terms µ 3 and µ 3 . We set the quartic couplings λ S and λ HS to zero since they need to be very small in any case to satisfy the XENON1T direct-detection bounds. We also set µ HS to zero as it does not contribute significantly to producing a first-order phase transition in the singlet direction in which we are interested. We also set the linear µ 3 1 term to zero to reduce unnecessary degeneracy, since it can be eliminated in favour of the cubic couplings, cf. appendix A.1. The parameter ranges we consider are m χ ∈ [30, 500] GeV, | sin θ| ∈ [0.01, 0.2], m 2 ∈ [10, 1500] GeV, |µ 3 |, |µ 3 | ∈ [0, 250] GeV and µ HS = λ S = λ HS = 0. To keep the scan denser, we employ a conservative perturbativity bound of π/2 on the absolute values of the quartic couplings.
To implement the constraints from the relic density, we use the micrOMEGAs code [69] with model files generated by the FeynRules package [70][71][72]. To search for the cosmological phase transitions and compute the nucleation temperature, the tunnelling action between the vacua, as well as the phase-transition quantities α and β, we employ the CosmoTransitions [73] code. These results were also checked by our own Mathematica code and using the FindBounce package [74].
In addition, we perform a separate scan to investigate the effect of the linear term µ 3 1 . In this case, we set µ 3 = µ 3 = µ HS = 0, so they are generated solely from elimination of the tadpoles. This elimination is achieved by shifting the field s → s + σ and then demanding that the linear term of the resulting potential vanishes. The effect of the shift on scalar potential parameters is given in appendix A.1. We parametrise this scan in terms of σ, not µ 3 1 to avoid solving a cubic equation. In this scan we consider parameter ranges |σ| ∈ [0.1, 200] GeV, |λ S | ∈ [0, 0.001], |λ HS | ∈ [0, 0.01]. We consider the singlet vev in the range |w| ∈ [1, 1.5 × 10 5 ] GeV and use it to fit the relic density to the observed value of Ω c h 2 within three standard deviations [43]. However, this linear term scan does not yield any points with strong first-order phase transition, nor is it particularly distinguishable in the direct-detection plots. For that reason we do not further discuss it separately.
The results of the scan in the plane of direct-detection cross section vs. DM mass are shown in figure 6. The points shown satisfy all theoretical and experimental constraints discussed in section 2 except for the direct-detection bounds by XENON1T which rule out the orange shaded region. The black and orange points produce a first-order phase transition while the gray points fail to do so. The black points are allowed by the XENON1T Figure 6: The direct-detection cross section vs. DM mass where the orange shaded region is excluded by the XENON1T experiment. The black and orange points produce a firstorder phase transition while the gray points fail to do so. The black points are allowed by the XENON1T bound while the orange points are excluded.
bound while the orange points are excluded. As is evident from the figure, there are ample regions of the parameter space which provide a viable DM candidate and lead to a firstorder phase transition. For these sets of parameters, we then determine the magnitude of the GW signal.
For the stochastic GW background signal, we recast the parameter points into the peak frequency-peak energy density plane, fixing κ t = 0.1. For the PIS curves, we assume the observational time t obs = 4 yr and threshold SNR ρ thr = 10 [16,62] for all the experiments. The most sensitive channel turns out to be bubble collisions for which we show the parameter points in the GW signal region in the left panel of figure 7. The black (orange) points are allowed (excluded) by XENON1T constraints.
We observe that obtaining parameter sets which lead to strong enough GW signal observable by future experiments becomes difficult in the generic parameter scan described in the beginning of this section. This is typical for a multi-dimensional parameter space constrained by multiple observables: the viable parameter space becomes more concentrated on lower-dimensional hypersurfaces and refined parameter scanning is needed.
Here, to explore the parameter space for the first-order phase transition more closely, we choose a few viable benchmark points from the general scan and search for more points in their vicinity along f rel = 1 hypersurfaces varying only two parameters at a time. This search strategy is illustrated in the right panel of figure 7, where we show the variation around a benchmark point   sets which correspond to models with sufficiently strong first-order transition to be visible in future searches for GWs. To better show the effect of the refined search over the full set of generated points, in figure 8 we show the points from the general scan discussed earlier in gray and the points generated by refined scanning in cyan.
To complete this discussion, we show how the scanned points are distributed by looking at various projections in the parameter space of the model. In figure 9 the gray points in both panels show all scanned points, while the orange (black) points show the points leading  to a first-order phase transition but are excluded (allowed) by XENON1T constraints on the direct-detection cross section. The cyan points are the ones generated by the refined scan.
In the left panel we show the dependence of the scanned points on the portal coupling λ HS and on the combination √ λ H λ S . The linear envelopes arise from the stability of the potential condition, while the curve bounding the points from above is due to the upper limit λ S = π/2 set by hand to guarantee a conservative bound on perturbativity and unitarity. As expected, the refined points cover very specific regions in the parameter space. Note that in the left panel some of the newly generated points go above the enveloping curve of the general scan. This is merely due to releasing the constraint λ S ≤ π/2 slightly but without endangering unitarity. The right panel shows the scanned points with respect to µ 2 S and µ 2 S . In figure 10 we show only the points which lead to a first-order phase transition. The orange (black) points are excluded (allowed) by the XENON1T experiment, while all other theoretical and experimental constraints discussed in section 2 are satisfied. The cyan points again correspond to the refined scan. The points are projected into the plane of the quantities relevant for the GWs signal. From the plots we see that both β and 1 − T n /T c quantifying the amount of supercooling are correlated with α, but the value of β and the value of the nucleation temperature are uncorrelated. The refined points follow the same correlation pattern, but the new points are more concentrated towards the region of large α.
In figure 11 the same points as in figure 10 are shown, but illustrating the dependence of GW signal on the parameters of the potential. In the left panel we see that the non-zero value of µ 3 allows for larger values of α. This is expected since a non-zero µ 3 contributes to a stronger phase transition. Similarly in the right panel we see that larger nucleation temperature requires a larger absolute value for the vev of the singlet field, |w|.
This concludes our analysis: we have established the parameter space of the model which provides for the observed DM abundance and is compatible with all present constraints from collider searches and direct and indirect DM detection. Moreover, this same parameter space allows for a first-order phase transition in the early universe whose resulting GW signal may be discovered in future observations.

Conclusions
In this paper, we have considered the most general model of pseudo-Goldstone DM that arises from the complex-singlet extension of the SM. Since the global U(1) symmetry is completely broken, the only remaining discrete symmetry is S → S * which stabilises the imaginary part of the singlet as a DM candidate.
In the Z 2 -symmetric case, in which the only U(1)-breaking term present is the mass term µ 2 S , the tree-level direct-detection cross section vanishes in the t → 0 limit. In the general case, we show that all other U(1)-breaking terms give a non-zero contribution to the direct-detection cross section in the t → 0 limit significantly increasing the interaction rates relevant for direct-detection experiments and lifting the protection due to the pseudo-Goldstone properties of the DM candidate.
However, we discovered that the symmetry-breaking parameters appear in certain combinations both in the pseudo-Goldstone mass and the tree-level direct-detection cross section. Setting such combinations to zero leads to cancellations which restore the pseudo-Goldstone properties of the DM candidate and suppress its direct-detection cross section in the t → 0 limit. Although running of the couplings upset this cancellation, the loop-level contributions can be mild and keep the direct-detection cross section moderately suppressed and still interesting from a phenomenological point of view. For these purposes we treat the model as a low energy effective theory where all symmetry breaking interactions may arise and are constrained by experiments and observations: in particular by the collider experiments and DM direct-detection experiments. We also consider constraints from indirect detection, but find that they do not presently constrain the parameter space significantly.
We also calculate the finite-temperature effective potential of the model and study the implications for phase transitions in the early universe. We consider a scenario where there is first a first-order transition in the singlet direction (w, v) = (w 0 , 0) → (w 1 , 0) with zero Higgs vev v, from which a second-order transition brings the fields to the electroweak minimum where both the singlet and the Higgs have non-zero vevs. The barrier in the first transition is generated by the singlet cubic couplings.
Although in most of the parameter space regions the first-order phase transition is not strong enough, we demonstrate that there are regions of the parameter space where a sizeable cubic coupling for a singlet can yield a stochastic GW signal with peak frequency of 10 −4 to 10 −2 Hz which may be detectable by future satellites BBO and DECIGO. Such couplings also tend to increase the direct-detection signal, which may be observable in future detectors such as the XENONnT experiment. The combination of cubics which suppresses the direct-detection cross section does not yield a strong enough GW signal.
couplings are shifted as A.2 Potential in the shift-invariant notation In ref. [75], the potential of the Higgs boson, h, and a real singlet, s, is given in terms of shift-invariant quantities. We can extend their formalism to write the potential in eq. (2.1) of h, s and χ in a shift-invariant way as where λ HR and λ HI are defined in eq. (2.7) and m 2 h , m 2 s and m 2 sh are the elements of the CP-even scalar mass matrix. A cosmological constant term appearing with a shift has been omitted.

B Annihilation cross sections
Here we define the tree-level annihilation cross sections of the pseudo-Goldstone DM, χ. We use short-hand notations s x ≡ sin x, and c x ≡ cos x for simplicity. It is also useful to define: and v rel = 2β χ . The annihilation cross section to fermionic final states is where the N c is 3 for quarks and 1 for leptons.

C Direct-detection cross section at one loop
In specific parts of the parameter space, the contributions from the symmetry-breaking interactions are suppressed at tree level. In such case the effect of loop corrections becomes relevant. Therefore, we will briefly discuss how to extend the analysis of ref. [28] to the case of cubic terms. Let us first briefly summarize the analysis of ref. [28] in the U(1)-invariant model with the Z 2 -symmetric mass term. We write the one-loop contributions to the direct-detection cross section at t → 0 limit as In the absence of other symmetry-breaking operators than the mass term for χ, λ 1L SI can be written as with short-hand notations s x ≡ sin x and c x ≡ cos x, and A 1 ≡4(m 2 1 s 2 θ + m 2 2 c 2 θ )(2m 2 1 vs 2 θ + 2m 2 2 vc 2 θ − m 2 1 ws 2θ + m 2 2 ws 2θ ), A 2 ≡ − 2m 4 1 s θ (m 2 1 + 5m 2 2 )wc θ − (m 2 1 − m 2 2 )(wc 3θ + 4vs 3 θ ) , A 3 ≡2m 4 2 c θ (5m 2 1 + m 2 2 )ws θ − (m 2 1 − m 2 2 )(ws 3θ + 4vc 3 θ ) . where A 6 ≡ 6 √ 2µ 3 vw + O(δ). (C.5) The functions A 0 , B 0 , C 2 and D 3 are the standard Passarino-Veltman functions; our definition of these agrees with the ones given in ref. [76]. The last term in eq. (C.4), O(δ), signifies that for δ = 0 several additional loop functions appear that cancel out in δ = 0 limit. Note that while the divergent parts of the loop functions of coefficients A 4 and A 5 cancel, the divergent part of the last term A 6 A 0 (m 2 χ ) does not. (Notice, however, that the divergent part vanishes in the limit m 2 χ → 0.) A new counter-term of the form 1 2 δµ 3 sχ 2 is needed, where (C.6)