Dark matter from a complex scalar singlet: The role of dark CP and other discrete symmetries

We study the case of a pseudo-scalar dark matter candidate which emerges from a complex scalar singlet, charged under a global U(1) symmetry, which is broken both explicitly and spontaneously. The pseudo-scalar is naturally stabilized by the presence of a remnant discrete symmetry: dark CP. We study and compare the phenomenology of several simplified models with only one explicit symmetry breaking term. We find that several regions of the parameter space are able to reproduce the observed dark matter abundance while respecting direct detection and invisible Higgs decay limits: in the resonances of the two scalars, featuring the known as forbidden or secluded dark matter, and through non-resonant Higgs-mediated annihilations. In some cases, combining different measurements would allow one to distinguish the breaking pattern of the symmetry. Moreover, this setup admits a light DM candidate at the sub-GeV scale. We also discuss the situation where more than one symmetry breaking term is present. In that case, the dark CP symmetry may be spontaneously broken, thus spoiling the stability of the dark matter candidate. Requiring that this does not happen imposes a constraint on the allowed parameter space. Finally, we consider an effective field theory approach valid in the pseudo-Nambu-Goldstone boson limit and when the U(1) breaking scale is much larger than the electroweak scale.


Introduction
In the last hundred years very strong evidence for the presence of non-baryonic cold dark matter (DM) in the Universe has been established [1, 2] (see also Refs. [3][4][5][6] for reviews). However, the Standard Model (SM) does not provide a suitable particle that could account for it. One of the simplest extensions of the SM that contains a DM candidate is obtained by adding a real scalar singlet, φ, equipped with a discrete symmetry that prevents its decay [7,8] 1 . The scalar potential of the SM is where H is the usual Higgs doublet, which after electroweak spontaneous symmetry breaking can be written in the unitary gauge as Here v = 246 GeV is the electroweak vacuum expectation value and h is the Higgs boson (in the absence of mixing with other scalars). This potential is now extended by adding three new couplings, where terms odd in φ are removed by imposing the discrete symmetry φ → −φ, which makes the φ particles stable if the symmetry is not broken spontaneously. The Higgs portal term, λ Hφ |H| 2 φ 2 , controls the annihilation of pairs of φ particles into SM states. The observed relic abundance can be reproduced via freeze-out for particular values of the DM mass (in the Higgs boson resonance, or for large values m φ 500 GeV [9][10][11][12]). The simplest model has been extensively studied [9,[12][13][14][15][16][17][18][19], and has been embedded in larger schemes that solve other problems of the SM such as neutrino masses, inflation and baryogenesis [20,21]. The mechanism is quite predictive because exactly the same coupling with the Higgs field that allows for DM annihilations gives rise to DM scatterings in direct detection (DD) experiments [22] (see also [23][24][25][26]) and, if the DM is light enough, induces invisible Higgs decays into DM [27][28][29]. Neither of the two effects has been observed and this sets important constraints on the model; in particular, DD experiments have left little parameter space for this scenario (see for instance Refs. [10,11,30,31]).
The potential above, which contains 19 real parameters, does not have any additional symmetry beyond the SM ones. However, the kinetic terms of the two new scalars do have an O(2) rotation symmetry and a shift symmetry φ 1,2 → φ 1,2 + c 1,2 . These symmetries can be used to remove the quadratic mixing m 2 12 φ 1 φ 2 in V II and the linear terms in V I (or the terms |H| 2 (ω 1 φ 1 + ω 2 φ 2 )). Therefore, the general potential adds 16 real physical parameters to the 2 parameters of the SM potential. Notice that these simplifications can always be used when the potential does not have any symmetry. However, if one imposes some symmetry that relates the parameters of the potential one should start from the most general potential compatible with the symmetry, check if there are parameters that can be removed by using the symmetries of the kinetic terms and then check whether the symmetries are broken or not by the global minimum of the potential. For instance, if we impose φ 1 ↔ φ 2 , we have that m 2 1 = m 2 2 , but m 2 12 = 0 is allowed; one could go to a different basis in which m 2 12 = 0 but with m 2 1 = m 2 2 , and in this basis the symmetry φ 1 ↔ φ 2 is hidden. Similarly, in the symmetry basis linear terms would be allowed, as long as δ 1 = δ 2 , and then one could ask about the conditions for this symmetry to be preserved by the vacuum.
If the potential does not have any symmetry the new particles do not carry any conserved quantum number that prevent them from decaying and, in general, they cannot be DM. The simplest way to have a sufficiently long-lived particle and therefore a good DM candidate is that it has at least one preserved symmetry. There are many symmetries one can impose; for instance, one can require that the potential is invariant under φ 1 → −φ 1 and φ 2 → −φ 2 . This will remove completely V III and the linear terms in the singlets. Then, if the symmetry is not spontaneously broken, the lightest of the two new particles will be completely stable. This model has been recently studied in detail in Ref. [32], where it has been shown that some of the problems of the model with just one scalar can be alleviated. One can impose many other symmetries: φ 1 → φ 1 , φ 2 → −φ 2 or φ 1 ↔ φ 2 or complete invariance under the symmetry of the kinetic term, O (2). Many of these symmetries lead to equivalent Lagrangians and some may or may not be broken spontaneously. It is then interesting to classify the different inequivalent possibilities that could lead to a good DM candidate and study the allowed parameter space for each of them. Since O(2) is isomorphic to U (1) times a reflection, it is more natural to explore the symmetries using the complex parametrization for the two real fields S ≡ (φ 1 + iφ 2 ) / √ 2. The possible models that could lead to a stable particle are quite obvious in this parametrization and can be classified in two groups depending on whether the complex field S acquires a vacuum expectation value (VEV) or not. In this work, we focus on the former case because it helps to avoid DD constraints [37][38][39][40][41]. Moreover it can also help to alleviate the instability problems of the scalar potential of the SM (see for instance Ref. [42]).
The paper is structured as follows. In Sec. 2 we discuss the general Lagrangian with two real scalar fields by using the complex parametrization, and we classify all the symmetries one can impose on it that could lead to a DM candidate. In Sec. 3 we study the scalar spectrum and couplings in the most general scenario. We outline the possible explicit symmetry breaking terms, which give rise to what we term as minimal models. In Sec. 4 we analyze the DM phenomenology, discussing the relic abundance, DD and Higgs invisible decays, presenting numerical results for a prototypical minimal model (the quadratic model) and comparing with the other ones. In addition, we discuss the possibility of a light DM candidate at the sub-GeV scale. In Sec. 5 we study the case with several symmetry breaking terms, which in general may lead to the spontaneous breaking of the symmetry that stabilizes the DM in large regions of the parameter space. In Sec. 6 we discuss an effective Lagrangian approach which is appropriate when the DM candidate is a pseudo-Nambu-Goldstone boson (PNGB). Finally in Sec. 7 we outline the main results of this work. We also provide several appendices with more technical details.

The complex scalar singlet Lagrangian and its symmetries
We write the complex scalar singlet field S as Then, the most general Lagrangian of the SM extended by S can be written as where L SM includes the SM potential, V SM , in Eq. (1), and V (H, S) contains all the interactions of the new complex scalar S. In the complex parametrization, the symmetry of the kinetic term is S → e iα S and S → S * , namely, U (1) times a reflection, which can be identified with a "dark" charge conjugation, or "dark" CP (DCP) since the kinetic term of the scalar also preserves parity. If the U (1) global symmetry is preserved by the Lagrangian but is spontaneously broken by the vacuum, there will be an exactly massless Goldstone boson and, of course, no DM candidate (it would contribute, though, to dark radiation). Therefore, if there is spontaneous symmetry breaking (SSB), in order to have a DM candidate we need both DCP conservation and, at least, one term that breaks explicitly the U (1). In this work, we are interested in this scenario. If the explicit symmetry breaking is much smaller than the VEV of S, the lightest of the scalars is a PNGB and can be the DM candidate. This limit is studied in Sec. 6.
In order to have a DM candidate the potential V (H, S) has to preserve, at least, a discrete symmetry that stabilizes the DM candidate. The possible discrete symmetries one can impose are the discrete subgroups of the kinetic term symmetry group, namely, O(2) ∼ U (1) + reflection, compatible with a polynomial Lagrangian. These are Z 2 (S → −S), Z 3 (S → e i2π/3 S) , Z 4 (S → iS) and the reflection, DCP, (S → S * ). 2 One may think that there could also be other discrete symmetries like S → −S * (φ 1 → −φ 1 , φ 2 → φ 2 in the real parametrization) or S → iS * (φ 1 ↔ φ 2 in the real parametrization), but all the symmetries that involve a reflection are physically equivalent to DCP, as can easily be seen by rephasing the field and redefining the parameters in the Lagrangian. Therefore, we define DCP as S → S * . A necessary condition for DCP conservation is that all couplings in the potential are real, like in the Standard Model CP. Notice, however, that since the potential only depends on |H| 2 this DCP can be defined independently of the standard CP. Then, it will not be affected by the SM violations of CP and will not induce additional CP violation in the SM.
Given the possible symmetries, it is natural to split the potential in the terms that preserve the U (1) symmetry, V 0 , and the terms that break it explicitly, V 1,full , V Z 2 ,full , V Z 3 and V Z 4 , classified according to the discrete symmetries they preserve. That is 2 If S does not get a VEV, all the symmetries discussed can be used to stabilize the DM. For instance, the case in which only Z2 is imposed, S → −S and general complex couplings are allowed has been discussed in Ref. [32]. In Ref. [43] the authors have explored the possibility of having asymmetric DM in an extended scenario with a real scalar singlet in addition to the complex scalar singlet. and The second, third and fourth terms respect a Z 2 , Z 3 and Z 4 symmetry, respectively. The parameters m 2 S , λ HS , λ S in V 0 are all real, therefore, V 0 is also invariant under S → S * . All the parameters of the symmetry breaking terms are, in principle, complex. Therefore, the general complex potential contains 3 real and 8 complex couplings, 19 real parameters, as in the real parametrization. The correspondence between the couplings in the real and complex parametrizations is given in Appendix A. If all the couplings in the complex parametrization are real, the 8 couplings m 2 12 , δ 2 , ω 2 , α 12 , µ 2 , µ 12 , β 12 , β 21 in the real parametrization vanish. Therefore, the presence of any of the latter implies that the DCP is explicitly violated in the potential.
By using a phase redefinition of S one can choose one of the couplings, for instance µ 2 S , real (and positive or negative), which would be equivalent to choosing m 2 12 = 0 in the real parametrization. Moreover, as in the real parametrization, the linear terms can be expressed in terms of other couplings when selecting the global minimum of the theory. Then, the general potential in the complex parametrization contains 4 real and 6 complex couplings, 16 real physical parameters in total, like in the real representation. From the 16 parameters, 3 preserve the U (1) symmetry and 13 break it. However, we should remark that if we are going to impose some symmetry, to count the number of physical parameters we should start with the general Lagrangian, including all the terms, and see how many of them survive the symmetry.
If the minimum of the potential is found for S = 0, U (1) undergoes SSB. Then, also Z 2 , Z 3 and Z 4 will be broken by the VEV of S. However, if all couplings are real, DCP will only be broken if the VEV of S is complex. It is easy to see that if there is only one U (1) symmetry breaking term its coupling can be taken real. Moreover, in that case, the VEV of S can always be chosen real [44], therefore DCP will be automatically preserved also by the vacuum. Alternatively, if there are several symmetry breaking terms, depending on the values of the couplings, the VEV can be complex even if the couplings are real and, therefore, the system can suffer from spontaneous DCP violation (see Sec. 5.1). Then, we arrive at the conclusion that if there is SSB of U (1) the only symmetry that could provide a DM candidate is DCP. In this case, the other discrete symmetries Z 2 , Z 3 , Z 4 can only be used to simplify the Lagrangian 3 and/or to motivate the presence of the DCP symmetry (e.g. in the case of Z 3 or Z 4 ).
The general case of SSB of U (1) with all possible explicit symmetry breaking terms contains many parameters and it is not very predictive. Moreover, the discrete symmetry DCP, necessary to have a DM candidate, must be imposed by hand, and even in this case it could be broken by the vacuum, leading to an unstable pseudo-scalar. On the other hand, if there is only one explicit U (1) symmetry breaking term, as discussed above, DCP is automatically conserved. Therefore, it makes sense to study first the cases with only one symmetry breaking term and then analyze how the addition of other terms modify the simplest scenarios. This can also be justified by using symmetries or by taking the softest symmetry breaking terms. We will make this analysis in Sec. 3.2.

The general case
Following the discussion above we will consider the general potential in Eqs. (11)(12)(13)(14)(15)(16) and require a DCP symmetry, e.g. that all the couplings in the potential are real. Then we will have 11 real parameters, 3 preserving U (1) and 8 breaking it. Notice that, when S is written in terms of real and imaginary parts, φ 1,2 , the Lagrangian is invariant under φ 2 → −φ 2 which is just the manifestation of the DCP symmetry. If unbroken by the vacuum, this symmetry will avoid φ 2 from decaying and it will be a good DM candidate. The real part φ 1 , however, has completely general couplings; in particular it has linear and cubic terms. Since the potential in terms of φ 1 is the most general potential we can write, it has the same form after a shift transformation φ 1 → φ 1 + σ which can be used to remove one of the parameters, for instance the cubic term µ H1 |H| 2 S or the linear term µ 3 S, leaving only 10 real parameters. However, for the time being, we maintain all of them, keeping in mind that if all of the DCP conserving couplings are present one of them is redundant. Then, we will assume that the scalar singlet also takes a VEV, v s . Using the linear parametrization for the singlet, we have The Higgs doublet VEV, v, can always be taken real, and to preserve the DCP symmetry we will require that also v s is real. Notice that the most general potential could lead to a complex v s , therefore breaking spontaneously the DCP and spoiling the stability of the DM candidate. In Sec. 5.1 we will comment on the conditions to avoid spontaneous DCP violation. Using the minimization equations for h and ρ , the bare mass parameters can be written in terms of the couplings and VEVs, Substituting them back in the potential allows us to compute the mass term of the fields. When all couplings and VEVs are real, θ does not mix with the other fields, and its mass is given by which displays the PNGB nature of the DM candidate θ, namely, its mass is zero if all the 8 symmetry breaking couplings vanish. On the other hand, the real parts of the fields, (h , ρ ), do mix with a mass matrix given by with the following matrix elements This mass matrix can be diagonalized by an orthogonal rotation of angle α where we have defined s α ≡ sin α and c α ≡ cos α in order to simplify the notation. The eigenstate h is chosen to be the 125 GeV boson observed at the LHC. In the numerical analysis, we take the mixing s α < 0.1 in order to satisfy experimental measurements of the Higgs signal strengths. One can trade-off some of the couplings in terms of the physical masses m h , m ρ and the mixing angle α Therefore, we can replace the 5 parameters of the symmetry preserving U (1) potential by Since the mass (squared) of the DM depends linearly on all symmetry breaking couplings one can replace one of them by the DM mass. In the simplified cases we will consider in the next section there is only one symmetry breaking term, therefore they will depend on 6 parameters, those in Eq. (29) plus the DM mass, m θ , from which four are unknown.

Minimal models
We introduce now minimal models, that is, with just one symmetry breaking term. In this case the coupling can be taken real, and therefore the DCP that stabilizes the pseudo-scalar is preserved.
The selection is motivated by either being the softest symmetry breaking term and/or by preserving a discrete symmetry. In V 1,full and V Z 2 ,full we just choose the softest term, that is, The potentials V Z 3 and V Z 4 in Eqs. (15) and (16) already contain only one term. At some point in the text we will refer to these minimal scenarios involving V 1 , V Z 2 , V Z 3 and V Z 4 as linear, quadratic, cubic and quartic models respectively.
In the numerical studies we impose the relevant theoretical constraints on the models: perturbativity, stability of the potential, and we have checked numerically that the minimum (v = 0, v s = 0) is the global minimum of the potential [37].
For the minimal models, we particularize the expressions for λ H , λ S and λ HS using the general Eqs. (25)(26)(27), with A = −1, 0, 1/3, 1/2 for the linear, quadratic, cubic and quartic models respectively. As can be seen in Eq. (32), the condition λ S > 0 is satisfied automatically in the minimal models except in the linear one, for which for the mixing angles we will consider, s α [10 −5 , 10 −1 ], it just reads m 2 θ m 2 ρ .

Radiative corrections
In principle the presence of a single explicit symmetry breaking term in the minimal models may induce radiatively some other breaking terms. For instance, let us consider the case of the quadratic model, with the potential in Eq. (31). At one loop it generates finite contributions to all couplings in V Z 2 ,full (Eq. (14)) and in V Z 4 (Eq. (16)). For m ρ m h , these corrections are of the order 4 Notice that λ H2 . This is also the case if the couplings come from higher-dimensional operators with spurions [47]. We have checked that the inclusion of the generated λ In the other minimal models (linear, cubic and quartic), the symmetry breaking terms do not generate any further ones. For the linear case this can be seen in the effective potential (i.e., it is not an interaction), while the cubic and quartic models in Eqs. (15) and (16) involve the only interactions allowed by a Z 3 and Z 4 symmetry, respectively, and therefore do not generate any other terms.

Dark matter relic abundance
In the Early Universe, the complex scalar singlet is in thermal equilibrium with the SM thanks to the interactions mediated by the Higgs portal coupling, λ HS . Then, at some point in the evolution of the Universe, the scalar S acquires a real VEV, breaking all possible symmetries except DCP, which, after SSB, is manifested by the pseudo-scalar θ as θ → −θ, which could therefore be a potential DM candidate. Once its interactions are slow enough compared to the expansion of the Universe, it freezes-out and thereafter its number density over entropy density, n/s, remains constant. Considering non-relativistic freeze-out, in large regions of the parameter space the interactions are too weak, the pseudo-scalar freezes-out too early, and the relic abundance is too large. However, there are parts of the parameter space that may produce large enough annihilations, so that the correct relic abundance is reproduced: 1. Resonances with the Higgs boson h or with the scalar ρ, which occur for 2m θ m h or 2m θ m ρ (See Fig. 1 (a)). In this region, kinetic equilibrium of the final states at freeze-out may not be a good assumption, and therefore there is some uncertainty in the parameter values for which the abundance is reproduced [48][49][50].
2. Direct annihilations into (lighter) pairs of scalars h and/or ρ, for m θ m h and/or m θ m ρ (See Fig. 1 (b)). The latter case is known as secluded dark matter (SDM). If λ HS = 0, for s α 10 −16 , ρ decays via mixing into SM states, as long as m ρ > 2m e .
3. Direct annihilations into somewhat slightly heavier pairs of hh, hρ, ρρ (See Fig. 1 (b)). This is known as forbidden dark matter (FDM). For m h m ρ m θ , annihilations into ρρ are normally larger and set the abundance. For m θ m h m ρ , which channel dominates depends on the mixing angle α. In Ref. [49] the case in which kinetic equilibrium of the final states at freeze-out may not be a good assumption has also been studied.

Non-resonant Higgs-mediated annihilations into SM states happening for DM masses above
100 GeV and at mixings s α larger than the previous cases. [51,52] (See Fig. 1 (a)).
Along this paper we use the DM relic abundance value from the Planck Collaboration [53] Ω obs h 2 = 0.120 ± 0.001 , where h is the dimensionless Hubble parameter, H = 100 h km/s/Mpc.

Direct detection
When there is non-zero mixing between the CP-even scalars, the pseudo-scalar may give rise to nuclear scatterings in underground detectors. In the limit of an exact Goldstone boson, the scatterings of θ are suppressed by the small momentum transfer. However, for the minimal models with different symmetry breaking terms, in principle significant interaction rates may be generated, which may allow one to distinguish them. We follow the analysis of Refs. [37,40]. The spinindependent DD cross section at tree level is given by where m N = 0.939 GeV is the nucleon mass, and f N = 0.28 (7) is the effective Higgs-nucleon coupling [54,55]. The effective DM-nucleon coupling λ SI reads as and the β i coefficients are defined by the interactions described in the Appendix B. One can see that in the zero momentum limit (t → 0), the effective DM-nucleon coupling goes as The expressions for the effective coupling in the small momentum limit are particularized for the minimal models in Table 1. As the coupling λ SI enters squared in the cross section, there will be no difference for the linear and cubic models, whereas in the case of the quartic model the effective coupling squared is a factor of 4 larger with respect to the former models. For the quadratic model the effective DM-nucleon coupling vanishes in the zero-momentum limit [37,40]. We will show this using an EFT approach in Sec. 6. We assume that the one-loop contributions are small compared to the tree level ones, as we are not in any particular point in parameter space where cancellations at tree level occur [37] (see also Refs. [56][57][58][59] for one-loop contributions in the case of the quadratic model, where indeed they are the dominant contribution).
In the numerical analysis, DD constraint from XENON1T [60] has been taken into account using the relative DM relic abundance of each point in the scan, and rescaling the experimental bound. We are therefore assuming that the local relic abundance scales like the global one, e.g.
ρ ∝ Ω. That is, every allowed point in the parameter space satisfies where σ XENON1T is the 90% confidence level upper limit on the DM-nucleon spin-independent cross section from XENON1T.

Higgs signal strength and invisible decays
For m θ < m h /2 and/or m ρ < m h /2 new decay channels of the Higgs boson into θ and ρ open up. The decay widths into these final states read with the β i coefficients defined in Appendix B. The Higgs invisible branching ratio is therefore given by where Γ SM h = 4.1 MeV is the SM Higgs decay width. The observed Higgs signal strength in the measured channels imposes a constraint on µ = c 2 α (1 − BR(h → inv)). A detailed analysis taking into account the contributions of both CP-even scalars to the different channels is beyond the scope of this work. However, in the limit of small mixing s α and very heavy scalar ρ, we can neglect the contribution of the latter, take c 2 α = 1 and directly impose the experimental constraints on the invisible Higgs width. If the mass of ρ is of the order of the Higgs boson mass or below (SDM), and if the mixing is very small, this is also a reasonable assumption. We therefore impose the 90% confidence level upper limit of BR(h → inv) < 0.16 obtained by the CMS Collaboration [61] (see also Ref. [62] for the ATLAS Collaboration results).

Numerical results
In the following we present the results of a numerical analysis of the scenarios discussed so far. First, we focus on the quadratic model as prototypical example, which also features suppressed DD rates, as we have discussed. The goal is to obtain all the regions where the correct relic abundance can be reproduced. Then we investigate the parameter space where it is possible to have a good DM candidate in each of the minimal models. Notice that if the symmetry breaking terms λ H2 and µ H1 in Eqs. (13) and (14) vanish, all the observables we consider depend on s 2 α . Therefore, we will only choose positive values for the mixing. Finally we analyze the possibility to distinguish the models comparing their predictions to different observables. The relic abundance has been computed using the code micrOMEGAs [63], see also Ref. [64].

The quadratic model
In this section we focus on one of the possibilities to obtain the correct relic abundance: setting the DM candidate θ in resonance with the scalar ρ, which corresponds to the option 1 discussed in Sec. 4.1. The relevant annihilation process is θθ → f f , where f refers to some SM particle. This setting allows one to expand the parameter space of DM masses, which are not necessarily needed to be close to the Higgs boson resonance [51,52]. We parametrize the deviation from the resonance condition with the dimensionless mass splitting parameter so that m ρ = (∆ + 1) m θ . In Fig. 2 we depict the (logarithm of the) mixing that reproduces the correct relic abundance in Eq. (34) as a function of the DM mass for different values of the dimensionless mass splitting ∆ and an illustrative value v s = 100 GeV. For masses below 4-5 GeV decays into hadrons are taken into account 5 and also constraints from having a light scalar mixing with the Higgs boson [65]. In the parameter space considered here, the most relevant constraints come from the invisible Higgs boson decays, the limits on rare B-meson decays, and being in thermal equilibrium with the SM particles in the Early Universe. The orange shaded region is excluded by the limits set on B → Kρ → K + "invisible" [66]. In our model, ρ decays into an invisible final state composed of two θ, and both escape the detectors leaving no signal. The calculation of the decay rate B → Kρ was performed using the expressions in Ref. [65]. The other relevant constraint is set by the invisible Higgs decay, shown as the blue shaded region. Notice that a detailed analysis of the relic abundance very close to the resonance has been performed in Ref. [48], and also for the case of a PNGB in Ref. [50], where significant differences appear in the computation of the relic density depending on whether kinetic equilibrium is maintained at freeze-out or not. Also, there is a huge sensitivity to the mass splitting, as can be seen by comparing the required mixing angles for the two different values of the parameter ∆. Therefore, Fig. 2 is intended to show that it is possible to have the proper relic abundance in the considered parameter space, but should be taken with a grain of salt regarding the precision of the exact value of the mixing angle for a given mass splitting that reproduces the abundance.

Comparison of minimal models
Now, we consider all the minimal models and study their allowed parameter space.
First, we focus on the restrictions that DD limits set on the minimal models, irrespectively of the relic abundance. In Fig. 3 we plot the allowed parameter space for the different models in the (m θ , m ρ ) plane for s α = 10 −1 for two values of v s (100 GeV in blue, 1 TeV in orange), after imposing the XENON1T limits [60]. In order to do that, we use the expression for the DD cross section in Eq. (35) and the λ SI values shown in Table 1. The quadratic model has a momentumsuppressed DD cross section at tree level and therefore it is not shown in the plot. Moreover, the effective coupling λ SI for the cubic model is similar to the linear case, apart from a relative sign, so we only plot the results for the latter one. For smaller values of the mixing, there are no restrictions from DD null-results, except for very light ρ masses. This can be understood by looking at the expressions of the effective coupling in Table 1 and its dependence on the ρ mass, which goes as λ SI 1/m 2 ρ for small m ρ values (but still larger than the typical momentum transfer at the XENON1T Experiment, O(MeV)).
In Appendix D we also discuss DD constraints for the symmetry breaking terms involving the Higgs field in Eqs. (13) and (14), which are qualitatively different. Now we impose the requirement of reproducing the relic abundance. First of all, we have checked that in all cases DM is in thermal equilibrium with SM particles in the Early Universe: in the resonance case the relevant process is θθ ↔ SM SM, and in the SDM/FDM scenarios the relevant ones are θθ ↔ ρρ, hρ, hh followed by ρ ↔ SM SM. Therefore, we can safely use micrOMEGAs [63] for the numerical computation. Moreover, in the case of SDM/FDM, we have checked that for the masses and mixings considered, the ρ particles always decay into SM states.
In Fig. 4 we show the results of a scan of the normalized mass splitting parameter ∆ and the mixing angle for three different values of the mass (for m θ = 40, 60, 130 GeV from left to right) and a fixed value of the VEV (v s = 100 GeV), where all the points satisfy that 0.5 Ω/Ω obs 1. As can be observed, the different freeze-out scenarios discussed in Sec. 4.1 are realized in separated regions: the Higgs resonance (h-res.) for m θ m h /2 at a fixed value of the mixing, the ρ resonance (ρ-res.) for ∆ 1, FDM when ∆ 0, SDM for ∆ < 0 and the non-resonant Higgs-mediated annihilations (non-res. h) for masses above 100 GeV and s α larger than the previous cases. It is important to point out that the Higgs resonance is only observed in the middle plot and that is in agreement with Fig. 2, where for m θ m h /2 the main DM annihilation channel into SM particles is mediated by the Higgs boson. As expected, the resonance regions are the ones changing the most among the different plots. Note that for the linear model, the proper relic abundance is reached only at the resonances. This is due to the fact that the theoretical constraint λ S > 0 restricts the values of the mass splitting to the region ∆ > 0. Therefore, SDM can not be realized in this case. Notice also that the perturbativity constraint (λ S < 4π) reduces the parameter space of each model in the case of large ρ masses. GeV from left to right, and v s = 100 GeV, for the minimal models described in the text. The green, blue, red and black colors correspond to the linear, quadratic, cubic and quartic models respectively. All these points fulfill the relic abundance condition 0.5 Ω/Ω obs 1. We also impose the XENON1T and the invisible Higgs decay constraints.
In Ref. [52] similar results where obtained for the case in which the relic abundance is reproduced through processes via resonances, and for the cases of large mixing and masses above 100 GeV. The former can be seen for all m θ masses considered in Fig. 4, and the latter in the case of m θ = 130 GeV, where the resonant region starts expanding from ∆ 1 (non-res. h). However the regions with proper relic values from FDM and SDM scenarios were not discussed in that work. In order to better understand the structure of strips of Fig. 4, we plot in Fig. 5 the relic abundance for fixed values of m θ = 60 GeV, v s = 100 GeV and s α = 10 −4 . We observe that for the case of FDM, close to ∆ 0, the linear model does not reach a cross section large enough in order to have the correct relic abundance. This is due to a partial cancellation happening among the diagrams (see the sign difference in Eq. (26) in the contribution to λ S of µ 3 with respect to the µ 3 or λ 4 one). In the quadratic model, for ∆ < 0 we see that SDM is allowed for two different values of the mass splitting. The explanation comes from the fact that the amplitude for the DM annihilations into ρ has two contributions, one proportional to m ρ and the other one to m θ , and in this model only the former is present, so when m ρ → 0 the amplitude vanishes at tree level 6 .
In Fig. 6 we display the mixing angle versus the DM mass setting the resonance mass condition in Eq. (41). We show plots for two different values of the VEV, v s = 100 GeV in the left panel,  Experimental constraints from invisible Higgs decays (blue), XENON1T experiment [60] (gray), the projection for XENONnT [67] (red dot-dashed line) and the thermalization condition (green) are shown. Same color code for the minimal models as in Fig. 4.
However, the opposite happens in the FDM and SDM scenarios, where the differences among models are clearly visible. This can also be observed in Fig. 7, where we plot the normalized mass splitting versus the VEV v s for the minimal models and different values of the DM mass (m θ = 40, 60, 130 GeV from left to right) with a fixed (small) mixing angle s α = 10 −5 . Notice that the SDM scenario can also happen for the cubic and quartic models when the value of v s is increased.  Fig. 8 we depict the results: the rescaled spin-independent DD cross section σ SI, resc = (Ω/Ω obs )σ SI for the different models, along with the current and future limits of the XENON1T experiment [60,67]. Notice that for the parameter values considered, in the case of SDM (∆ = −0.1, for instance) the DM is under-abundant, so we do not show this case. Remember that for the quadratic model in the zero-momentum limit there is an exact cancellation at tree level in the DD cross section (see Table 1). Therefore, its dominant contribution is at one loop and is very suppressed. We can see that in the FDM region (∆ = 0.1), assuming astrophysics are under control (say, the standard halo model), a precise-enough positive measurement of a DD signal could allow one to distinguish among the minimal models. This is not the case for the region close to the ρ resonance (∆ = 1.1), and also close to the Higgs resonance (m θ 60 GeV), where the required parameters are quite insensitive to the minimal model considered.  1 (right, ρ resonance). Constraints from perturbativity and invisible Higgs decay have been taken into account. We also plot in gray the exclusion region from the current XENON1T experimental limit [60] and its projection [67] as a red dot-dashed line.
We have checked that in all the minimal models there is a significant temperature dependence of the thermally averaged annihilation cross section in the resonance scenario, allowing to avoid indirect detection (ID) bounds [68]. This was also found in the case of the quadratic model in Ref. [51]. However, in the SDM regime the variation on the temperature is more subtle (see next section). Of course, ID bounds do not constrain FDM, given the absence of DM annihilations at zero temperature in this case.
Finally, we have also studied the DM self-interacting cross section for each model and found that the maximum values were reached in the resonance with the scalar ρ. However, in the parameter space considered the values of the cross section were much smaller than the ones constrained by clusters, 1 cm 2 /g [69].

Light dark matter
In this section we discuss the possibility of light DM (e.g. at the sub-GeV scale) in the minimal models 7 . Let us analyze in the following the different freeze-out scenarios described in Sec. 4.1. We impose that the DM candidates are in thermal equilibrium in the Early Universe. Furthermore, constraints from DD experiments, namely CRESST-III and DarkSide-50 [71,72], and invisible Higgs decays are applied.
1. In the ρ resonance, the lightest possible mass for the DM candidate can be read from Fig. 2 and it could be at the sub-GeV scale. As discussed in Sec. 4.4.2 (Fig. 6), it is difficult to disentangle the minimal models in the resonance condition, therefore the same lower bound for the DM mass applies for the rest of the models. Constraints on the mixing angle coming from meson decays into invisible states (B → Kρ and K → πρ [65]) forbid lower DM masses.
2. In the SDM/FDM scenarios, restrictions on the mixing angle come basically from existing limits on light scalars that mix with the Higgs boson (see Fig. 8 in Ref. [65]). For the SDM regime, ID bounds must also be considered.
As an illustrative example, in Fig. 9 we take m θ = 1 GeV and perform an scan for the minimal models where the points fulfill the condition 0.5 Ω/Ω obs 1: on the left panel we plot ∆ versus v s with s α = 10 −5 ; on the right panel we depict ∆ versus s α with v s = 10 GeV. For this value of v s , the constraint from the invisible Higgs decay almost excludes the ρ resonance solution, and ID bound from Fermi-LAT [68] forbids SDM for ∆ −0.1.

Beyond minimal models
Here we will try to go a bit beyond the minimal models discussed above, which involved just one dominant symmetry breaking term. We do so by combining pairs of minimal models, that is including simultaneously two symmetry breaking terms. The main goal is to study if the parameter space opens up. However, notice that the inclusion of the additional couplings (even if real) can lead to SSB of DCP and spoil the DM candidate, as it would not longer be stable. Therefore, we should require that this does not happen, which leads to important constraints on the parameter space of these non-minimal models.

About spontaneous dark CP violation
In this section we study the case in which the potential contains two explicit U (1) symmetry breaking terms, which we write as where n, n = 1, 2, 3, 4 correspond to the linear, quadratic, cubic and quartic model interactions, andã andb are generically the couplings of the two interactions we consider. Without loss of generality we take n < n . DCP conservation, defined as S → S * , requires that the two couplingsã andb are real. To avoid SSB of DCP we should also require that the VEV of S is invariant, S = S * , and therefore, S must be real. Moreover, since the rest of the Lagrangian is invariant under phase transformations, we can make a redefinition S → −S and take always S real and positive. Then, to explore the conditions under which DCP is not spontaneously broken it is convenient to use the exponential parametrization (see Refs. [44,73]) because in this parametrization the field G only appears in the symmetry breaking part of the Lagrangian, Eq. (42). Then, according to the discussion above, we require that the global minimum of the potential is found at σ = 0 (this will fix the value of v s in terms of all the parameters of the potential) and G = 0 so that √ 2 S = v s > 0. In this parametrization the potential contains the terms where we have defined the dimensionless couplings which are expressed in terms of the original couplings of the Lagrangian as given in Table 2.
n, n Minimal model a, b Table 2: Relation of the effective parameters a and b defined in Eqs. (42,(44)(45) and the explicit symmetry breaking couplings in the potential for the minimal models.
Then, to avoid spontaneous violation of DCP one has to require that G = 0 is the global minimum of V sb in Eq. (44). Obviously, the first derivative of V sb in G = 0 is zero, while the second derivative is −v 2 s (n 2 a + n 2 b). Therefore, to have a local minimum at G = 0 one should require b ≤ −a(n 2 /n 2 ). Then, one has to check that it is the global minimum of V sb (G) by comparing it with other minima, which has to be done case by case. We thus find that the conditions for DCP conservation are a ≤ 0 and These constraints are displayed in Fig. 10 in the (a, b) plane and, when expressed in terms of the explicit symmetry breaking parameters of the general potential in Eqs. (11)(12)(13)(14)(15)(16), must be added to the theoretical constraints mentioned in Sec. 3.2. As an illustrative example of the consequences of combining two minimal models, in Fig. 11 we plot the results in the planes (s α , ∆) on the left and (v s , ∆) on the right for having a suitable DM candidate when considering the cubic and quartic symmetry breaking terms, both alone and together. The results show that including both terms allows values of m ρ and v s in between the individual cases (minimal models); that is the region covered by the gray dots in the figure. Analogous results were found when considering the other combinations of minimal models.
In conclusion, having two symmetry breaking terms enlarge the parameter space for a good DM candidate to regions bounded by the minimal models. This is a consequence of the additional constraints required to avoid the spontaneous violation of DCP shown in Fig. 10.   The pseudo-Nambu-Goldstone Boson limit: an EFT approach 6

.1 Effective operators
For the following discussion it is useful to consider the exponential parametrization, Eq. (43), in terms of the radial mode σ and the angular mode G. Let us now assume that S takes a large VEV as compared to the explicit symmetry breaking terms. Then, the mass of the angular mode G is much smaller that the symmetry breaking scale m G v s and it can be considered a PNGB. At low energies, the scenario which involves a PNGB G stabilized by a discrete symmetry (DCP) G → −G can be parametrized by where the PNGB mass m G , the quartic coupling λ G and the Higgs portal λ HG break the shiftsymmetry. These correspond at high energies to linear, quadratic, cubic and quartic terms in the complex singlet scalar S in the linear parametrization. Notice that they are not suppressed by derivatives, e.g., they can give significant contributions for instance in DD. Let us parametrize the breaking with the charge n of the field S (e.g., V ∝ S n + H.c.). We can now proceed to integrate-out the radial mode σ . In terms of the UV parameters, the Wilson coefficients read and m G = m θ is given in Eq. (20). Notice that m G and λ G (and c G and λ HG ) are related with each other. We have expressed these coefficients in terms of physical parameters using also that, in the limit v s v, If one considers a certain shift symmetry G/v s → G/v s + 2π/n at low energies, the non-derivative terms may be dropped and the allowed terms in the potential take the form where and we have assumed a renormalizable UV completion, which sets the upper limit in the sums. These d n (with n = 1, . . . , 4) low-energy terms have mass-dimension 4 and correspond at high energies, respectively, to terms linear, quadratic, cubic and quartic in the complex singlet scalar S in the exponential parametrization, e.g., respecting respectively the discrete symmetry DCP, Z 2 , Similarly, the c n (with n = 1, 2) have mass dimension two and correspond at high energies, respectively, to terms linear and quadratic with the Higgs doublet, e.g., respecting respectively the discrete symmetry DCP and Z 2 . Once expanded, they yield the mass and quartic interactions of the PNGB. It is interesting to note that if we integrate by parts the derivative interaction in Eq. (47) we obtain where in the last step we use the Klein-Gordon equation of motion for G, which is correct for G on-shell, as in the case of DD experiments. If we substitute back in Eq. (47) we find an additional contribution to the Higgs portal coupling, which exactly cancels for n = 2 and produces a suppression in DD experiments, since the other terms always involve the momenta of the PNGB. This cancellation was already noticed in Ref. [37].
Regarding the relic abundance, there are two options depending on the scale of v s : • If v s v, unless very small couplings are involved, we expect m h m σ , m G . In order to explain the relic abundance, the DM G needs to be at the radial mode (σ ) resonance, e.g., v, or almost degenerate with the σ , e.g. m G m σ v (for SDM/FDM). In either case, the small mass splitting between σ and G is not naturally achieved, as the symmetry breaking terms need to be comparable to v s , and therefore, strictly speaking, in that case G cannot be considered a natural PNGB.
• If v s v, unless very small couplings are involved, we expect m σ m h , m G . We can integrate-out the radial mode σ , and, in this case, the prediction is that the PNGB needs to be at the Higgs resonance to explain the relic abundance, e.g., m G m h /2, see Refs. [40,51], or almost degenerate with the σ , e.g. m G m σ v (for SDM/FDM). In this scenario, the small mass splitting between h and G is achieved for breakings of the order of the electroweak scale (which is below v s ). We can use an effective Lagrangian with just the PNGB and the Higgs field (see also Refs. [74,75]), and compare to the results of the complete model (for instance the relic density and direct detection).

Numerical results
In Fig. 12 we plot the relic abundance versus the DM mass for both the effective Lagrangian of the PNGB in Eq. (47) and the full quadratic model which includes the radial mode, with the symmetry breaking term in Eq. (31) that gives the mass to the DM candidate.
In the complete model (blue dashed line) we have set s α = 0.1 and v s = m σ = 10 3 GeV, where m σ is the mass of the radial mode. In the effective Lagrangian approach (orange line), we use the Wilson coefficients c G and λ HG from Eqs. (48) and (49) with the values for s α and v s mentioned before. Notice that, as we are comparing the effective Lagrangian with the full quadratic model, n = 2 in the λ HG coefficient. We can see both the Higgs and the radial resonances in the case of the complete model. The EFT reproduces very well the full model, including the Higgs resonance, up to DM masses below v s /6, above which the EFT cannot be trusted. In the full model, we also see the solution at m G m σ 1 TeV (forbidden/secluded, annihilations into σ bosons), which of course cannot be captured by the EFT. In the full model and in the EFT we see a small kink at m G m h corresponding to the opening of the annihilation channel GG → hh. A small kink at m G = 1 TeV, corresponding to the opening of the channel GG → σ σ , can also be seen in the full model.  31), which also includes the radial mode σ (blue dashed line). The correct value for the relic abundance is shown as a black dotted line, e.g., the region above is excluded and in the region below the PNGB is under-abundant.

Conclusions
A real scalar singlet, stabilized by a discrete symmetry, is one of the simplest candidates for DM. However, by combining relic abundance constraints with direct detection null-results, the allowed parameter space has been almost completely ruled-out, leaving only a few allowed spots: at the Higgs boson resonance, or at high masses, where annihilations into Higgs boson pairs opens up.
In this work we have extended this minimal scenario in the simplest possible way: a complex scalar singlet, charged under a global U (1) symmetry. If the symmetry is not broken, neither spontaneously nor explicitly, the allowed parameter space is very similar to that of the real scalar singlet 8 . If the symmetry is preserved at the Lagrangian level but not respected by the vacuum of the theory there will be an exact Goldstone boson in the spectrum and, since it is massless, it is dark radiation and cannot constitute the DM of the Universe. Therefore, the symmetry must also be broken explicitly at the Lagrangian level, yielding a pseudo-Nambu-Goldstone Boson as the DM candidate. The explicit symmetry breaking of the U (1) needs to preserve at least the dark CP that stabilizes the DM particle. This is the case we have studied in this work.
We have considered 4 cases with only one explicit breaking term, which we call linear, quadratic, cubic and quartic. The choices are motivated by either being the softest possible symmetry breaking term, or by a discrete symmetry. All these models preserve a dark CP symmetry, while the quadratic, cubic and quartic models are also invariant under Z 2 , Z 3 and Z 4 , respectively. These are the simplest models in the sense that the breaking involves only the new complex scalar singlet and it is the softest within each class. Furthermore, they are quite stable under radiative corrections, as only the case of the quadratic model generates, at one loop, further symmetry breaking terms, which are, however, suppressed.
The models are necessarily very predictive, with just four new parameters: m θ , m ρ , s α and v s . However, the allowed parameter space to obtain the correct relic abundance is enlarged significantly with respect the case of only one real scalar because the presence of an additional scalar, ρ, that mixes with the Higgs boson. Then, the allowed annihilation channels are now duplicated: being at a resonance with h and/or ρ, or at high masses, where annihilations into h and/or ρ pairs may open up. The last case, when h and/or ρ are (a bit) more massive than the pseudo-scalar is known as forbidden DM. Moreover, the relic abundance can also be reproduced in the limit of small mixing, via annihilations into h and/or ρ pairs. Annihilations in ρ pairs when m ρ < m θ provide an example of secluded DM.
Our analysis shows that these minimal models may potentially be distinguished among themselves if a positive signal in direct detection is observed. Measuring the DM mass (say, by a gamma ray line in indirect detection), for instance, can yield further information regarding the underlying symmetry. Moreover a positive detection of a new scalar would yield further information on the mixing, s α , and on the symmetry breaking scale. It should have couplings to SM fermions like those of the Higgs boson but suppressed by the mixing with respect to the latter, with possible extra decay modes ρ → θθ and/or ρ → hh if kinematically allowed.
In addition to this, we have analyzed the possibility of light DM. At the ρ resonance the lowest DM mass could be at the sub-GeV scale where annihilations into muons and hadrons are still possible, whereas in the forbidden/secluded DM cases the mass could be even lower, specially in the forbidden scenario for which indirect detection bounds do not apply.
We also study the case in which two of these symmetry breaking terms appear simultaneously, to see if the parameter space increases. In this case, the dark CP symmetry, which for just one symmetry breaking term (that can always be taken as real) was preserved also after SSB of the U (1), may be violated by the vacuum. We have derived the restrictions for all possible pair combinations of the symmetry breaking terms that appear in the minimal models so that the dark CP is preserved after SSB and, therefore, the pseudo-scalar is stable. We find that, once these DM stability constraints are taken into account, the allowed parameter space opens up precisely in the region between the two minimal models. Therefore, we can conclude that adding more symmetry breaking terms (at least by pairs) fills up the parameter space between the minimal models.
Finally, for small explicit symmetry breaking compared to the scale of SSB of the U (1), we have obtained an effective low energy Lagrangian including both the usual Higgs portal as well as a derivative Higgs portal. This effective Lagrangian turns out to be very useful in order to analyze DM phenomenology, specifically the relic abundance in the case of the Higgs boson resonance and direct detection.

C Hadronic decay modes
We follow the procedure described by Ref. [12] with the cross section for the ρ-mediated 9 s-channel DM annihilations into hadronic final states written as .
The β ρθθ coefficient is given by the Lagrangian in Eq. (56), the decay width of the ρ into hadrons, Γ ρ→hadrons , is taken from Fig. 4 in Ref. [65], and the full ρ width, Γ ρ,full , is just the Γ ρ→θθ , which is the dominant one in the considered parameter space. This decay width is described by

D Direct detection for explicit symmetry breaking terms with the Higgs field
Here we also consider the following explicit symmetry breaking terms with the Higgs field, V H2 = 1 2 λ H2 |H| 2 S 2 + H.c. .
We focus on the DD constraints (irrespectively of the relic abundance, as in the first analysis in Sec. 4.4.2) for models including these potentials. We provide the effective DM-nucleon couplings generated in Table 5.  Table 5: Effective DM-nucleon coupling that enters in the DD cross section in terms of the physical parameters v, v s , m h , m θ , m ρ and s α . Fig. 13 shows the allowed parameter space by the XENON1T bound for both models. We plot two values of the VEV (v s = 10 2 , 10 3 GeV in blue and orange, respectively) and set s α = 10 −3 , which is the typical value for the mixing in order to get the correct relic abundance near the resonances (see Fig. 6). We can see that DM masses larger than 10 GeV are excluded by the experimental constraint from XENON1T.  Figure 13: Parameter space that is allowed by the XENON1T bound for the models with the symmetry breaking terms in the potential V H1 and V H2 in Eqs. (61) and (62). Blue (Orange) colored region corresponds to v s = 100 (1000) GeV. Black and gray dashed lines represent the resonance condition with the ρ (m ρ = 2m θ ) and the degenerate case (m ρ = m θ ) respectively.