Renormalization Group Study of the Minimal Majoronic Dark Radiation and Dark Matter Model

We study the 1-loop renormalization group equation running in the simplest singlet Majoron model constructed by us earlier to accommodate the dark radiation and dark matter content in the universe. A comprehensive numerical study was performed to explore the whole model parameter space. A smaller effective number of neutrinos $\triangle N_{eff}\sim 0.05$, or a Majoron decoupling temperature higher than the charm quark mass, is preferred. We found that a heavy scalar dark matter, $\rho$, of mass $1.5-4$ TeV is required by the stability of the scalar potential and an operational type-I see-saw mechanism for neutrino masses. A neutral scalar, $S$, of mass in the $10-100$ GeV range and its mixing with the standard model Higgs as large as $0.1$ is also predicted. The dominant decay modes are $S$ into $b\bar{b}$ and/or $\omega\omega$. A sensitive search will come from rare $Z$ decays via the chain $Z\rightarrow S+ f\bar{f}$, where $f$ is a Standard Model fermion, followed by $S$ into a pair of Majoron and/or b-quarks. The interesting consequences of dark matter bound state due to the sizable $S\rho \rho$-coupling are discussed as well. In particular, shower-like events with an apparent neutrino energy at $M_\rho$ could contribute to the observed effective neutrino flux in underground neutrino detectors such as IceCube.


I. INTRODUCTION
In two previous studies [1], [2] we have constructed extensions of the singlet Majoron model [3] [4] with the motivation of accommodating possible new relativistic degree of freedom commonly known as dark radiation (DR) in cosmological models and also to provide a viable dark matter (DM) candidate. The Majoron which is the Goldstone boson from the spontaneous breaking of the global U(1) ℓ lepton symmetry is identified as the DR. This breaking is facilitated by a Standard Model (SM) singlet carrying lepton number of two units. We then add a non-Higgs singlet complex scalar with lepton number ℓ = 1. After symmetry breaking a stable scalar DM is obtained. This model preserves the simplicity of the Majoron model and connects the Type I seesaw mechanism to the dark sector which consists of dark matter and perhaps dark radiation. Moreover, the main motivation was to study the physics consequences of identifying the Majoron as dark radiation. In particular if the decoupling temperature T dec is at the muon mass, m µ it will give a contribution to the effective relativistic degree of freedom ∆N eff = .39 which is a sweet spot pointed out in [5] yielding N eff = 3.44. This is higher but not inconsistent with the 2015 Planck result of N eff = 3.15 ± .23 at 1σ C.L. [6]. When other data such as South Pole Telescope and Atacama Cosmology Telescope results are included the central value is higher. Although the Planck 2015 data is consistent with the SM prediction, but the statistical significance to rule out DR is still very poor. For example, △N ef f = .33 is still consistent with the Planck data at 1σ C.L. Since a massless Majoron is automatically built in in this model which implements a spontaneously-broken global U(1) l and it always contributes to the dark radiation. The relevant question is to determine how much can it contribute to △N ef f . This prompted us to reexamine the Majoron dark radiation model by allowing T dec to be higher than m µ which will reduce ∆N eff . For example for T dec around 2 GeV ∆N eff = .05 since now more degree of freedom contributes to the energy density. A consequence pointed out in [1] and [5] when taking T dec at m µ gives rise to a scalar of mass less than a few GeV. We find that increasing T dec will raise the mass of this scalar to the tens of GeV range. This will in turn change the impact on Higgs boson decays since this scalar in general mixes with the SM Higgs boson.
In a recent study [7] the impact of vacuum stability on the singlet Majoron model with high scale Type I seesaw model for neutrino masses was investigated. Besides the high seesaw scale the mechanism also requires Yukawa couplings of the righthanded Majorana neutrinos to active neutrinos via the Higgs field. Spontaneous electroweak symmetry breaking gives rise to a Dirac mass term. This is the second crucial ingredient of type I seesaw mechanism.
Hence, for it to work the electroweak vacuum must be stable for when lepton number breaking occurs. In other words the seesaw scale must be lower than the energy, µ SM V S , where the electroweak vacuum becomes unstable which is known to be around µ SM V S ≃ 10 10 − 10 12 GeV [8][9] [10]. For the Majoron model lepton number is spontaneously broken and hence the stability of the singlet scalar that breaks this symmetry must also be taken into account.
It was found that the stability of the SM can be extended to the GUT scale without invoking metastability. See [11] for a similar discussion on the vacuum stability by identifying the axion as DR.
In this paper we study how RG considerations impact the parameters of dark matter and dark radiation sector of the Majoron model. We calculate the one loop beta function for the renormalization group running of the all the relevant parameters of the Majoron dark radiation model of [1]. We find that the stability of the scalar vacua has very important effect on the parameters of the theory. In particular the dark matter candidate ρ will have a mass in the range of the lepton number violation scale; i.e. in the several TeV range. This is vastly different from the usual studies which did not take into account scalar vacua stability.
We find that in doing so can lead to interesting astroparticle physics consequences. Since the DM is heavy and there is a light scalar in the spectrum, bound state of DM can be formed if the triple scalar coupling is strong enough. We show that this can indeed take place in a large region of the model parameter space. We speculate that DM annihilation into two Goldstone bosons will be enhanced. One would then expect a Goldstone component in the high energy cosmic ray spectrum.
We organize the paper as follows. In section two we give a summary of the model and the RGEs of the relevant parameters. It is sufficient to use the 1-loop result for the beyond SM physics. This is followed by details of the numerical study of the model including the solutions of the RGEs. In section IV we discuss the phenomenological consequences of the results we obtained. Finally we conclude in section V.

II. THE MODEL
A singlet Higgs field S which carries lepton number ℓ = 2 and a non-Higgssed scalar field Φ with ℓ = 1 are added to the particle contents of the SM. Realistic implementation of the Type-I seesaw mechanism will require adding at least two singlet Majorana righthanded neutrinos N Ri , i = 1, 2, to the SM. Since the details of the neutrino physics such as active neutrino masses and oscillations are not relevant to this study and we can just take one righthanded neutrino for simplicity without affecting the physics we are interested in. Extending to the realistic case of two or more righthanded neutrinos is straightforward.
Also,the SM Higgs field is denoted by H.
Due to the U(1) l symmetry, Φ will not have a trilinear coupling with H;thus, it will not contribute to the Majorana mass of N R . Its Dirac mass type of couplings to the active neutrinos are also forbidden since it is an SU(2) singlet. Therefore, much of the Majoron model is not changed and its simplicity is retained.
The most general scalar lagrangian is given as and we take κ to be real and m 2 Φ > 0 so that Φ = 0. Using the usual linear representation of scalar fields we expand them as follows and use the U-gauge for the Higgs The physical fields areŜ = (h, s, ρ, χ) and ω is the massless Goldstone boson named the Majoron. With this one obtains the scalar mass matrix squared: Note that κ splits the degeneracy of the ρ and χ masses and we require . We take ρ to be the DM and its stability is guaranteed by Z 2 dark parity which remains after spontaneous symmetry breaking of U(1) ℓ [1].
In terms of component fields the scalar potential becomes whereκ ≡ λ ΦS v s + κ.
It is clear that (h, s) are not yet mass eigenstates denoted by (h 1 , h 2 ). They are related by the usual rotation: with the mixing angle θ given by We shall identify h 1 ≡ H as the SM Higgs which has a mass of 125 GeV. Note that in the small mixing limit, and h 2 ≈ S. The Lagrangian responsible for the seesaw mechanism is given by where L = (n L , e L ) T is the SM lepton doublet andH = iσ 2 H * . After symmetry breaking we get We wish to identify ω as the DR and the amount it contributes to ∆N eff depends on when it decouples from the thermal bath. In particular we are interested at decoupling temperatures around QCD phase transition, charm mass and tau mass; then ∆N eff = 0.055, 0.0451, 0.0423 respectively. The effective Lagrangian for ωω → ff is given by where f denotes the SM fermion in the thermal bath. The rate of ff ↔ ωω is estimated to where N f c is the color of f . Apparently the specie, which will be just denoted as f , with the largest values of (m 2 In order for ω to play the role of DR the collision rate of ω into a pair of fermions must be approximately the Hubble expansion rate at T dec , thus; For T dec at around m τ , both charm and tau must be considered and we have N c m 2 f → N c m 2 c + m 2 τ ; otherwise m f is the mass of fermion nearest to the T dec . In Eq.(12) the parameter λ HS is controlled by the mixing of the scalar S with the Higgs. It is expected to be small. Eq. (12) shows that M s is in the 100 GeV range if T dec is around the charm mass. This is to be compared to decoupling at m µ which leads to a light scalar of mass less than a few GeV. Details for the latter case can be found in [1].
For a given T dec the largest M S can be estimated by considering the Higgs invisible decay width, which is experimentally given by Γ inv H < 0.8 MeV [12,13]. Since the Majoron is massless, the SM Higgs can always decay into two Majorons and this width is denoted by Γ ωω . There may be other invisible modes available; thus Γ ωω ≤ Γ inv H . Moreover, Γ ωω is given by in terms of the physical mass of the SM Higgs. Using the relations between sin θ and v S and the mass eigenvalues plus the decoupling condition, Eq.(12), we can rewrite the above as From the above inequality, the upper bond of M S can be easily solved analytically. We should just denote the solution as M max S (T dec ), shown in Fig.1(a), since the explicit form is not important. A direct search for the light neutral scalar denoted by S here, at OPAL [14] yields an upper limit of the size of mixing between H and S. On the other hand, the mixing will modify the SM Higgs coupling to anything by a cos θ factor. At the LHC, the signal strength µ f i for a specific production and decay channel i → H → f is defined as If the SM Higgs invisible decay width Γ inv H ≪ Γ SM H , which is the case in our numerical study, then BR f ≃ (BR f ) SM and µ f i ≃ cos 2 θ is predicted in our model. The best-fit of signal strength µ = 1.1 ± 0.11 is given in a recent ATLAS and CMS combined global analysis on all production process and decay channels with data taken at √ s = 7 and 8 TeV [15]. This indirect bound amounts to sin 2 θ < 0.13 at 2 σ level which has also been implemented for M S > 60 GeV in our numerical study, Fig.1(b). This latest bound is derived from LHC Run-1 data only, and the expected sensitivity of Run-2 will be discussed in the phenomenology section later.
Although many of the parameters in the scalar potential Eq.(5) are unknown, we can gain some information by demanding that the stability of the scalar sector in the appropriate range. The stability of the electroweak vacuum govern by the sign of λ is well known to be at best metastable [8,9] for the SM. On the other hand singlet scalars generally helps to stabilize the electroweak vacuum. In our model we also require that both S and Φ should have stable potentials for consistency reasons. The scales of stability are given by the RG running of the parameters. The RGEs for the SM couplings are easily found in the literature and we will not repeat them. The relevant RGEs for the new parameters calculated with 1-loop β functions are given below: where y t , g 2 , g 1 are the t-quark Yukawa coupling, SU(2) and U(1) Y gauge couplings respectively and t ≡ ln Q 2 with Q 0 an arbitrary renormalization point. We have omitted all light fermion Yukawa couplings including those of the active neutrinos since they are all very small. The running of the y ν 's can be shown to be unimportant for us [7]. For stability we require λ i > 0 and λ ij > −2 λ i λ j where i, j = H, S, Φ. The RGE's for these couplings by themselves are not sufficient to determine whether any of them will turn negative at high enough energies. We need boundary conditions at some lower energies. Since M S < 100 GeV, we choose this scale to be m Z . The values of the couplings at this scale will be given by numerical scan that has to satisfy other constraints we impose on the model. Details are given in the next section.
One important input comes from DM considerations. Our addition to the minimal Majoron model produces a WIMP DM candidate. Due to a Z 2 dark parity the lighter of ρ and χ will be the DM. Without loss of generality we take that to be ρ and their masses are split where for i = W, Z, H, f, S, χ, and the subscripts denote the final state. The couplings in the scalar mass eigenstates are given as At high temperatures these will give σv and for the correct relic abundance the total σv should be ∼ 3 × 10 −26 cm 3 /s. A numerical scan is performed as described in the next section to obtain possible values of unknown couplings.
Next we discuss the constraint imposed by the limits from direct DM detection since there is no convincing signals yet. In our model the scattering of ρ off the nucleus of the detector will deposit energy. The scattering ρ + n → ρ + n where n denotes a nucleon proceeds via the t-channel exchange of H and S. It is often to parameterize the SM Higgsnucleon-nucleon coupling by ηg 2 M n /(2M W ) [16] where M n is the nucleon mass and η is a parameter represents the uncertainty in the coupling. In the interaction basis the h n n and s n n couplings become c θ ηg 2 M n /(2M W ) and s θ ηg 2 M n /(2M W ) respectively. And the tree-level cross section in terms of physical masses of H and S is where the reduced mass is We take η = .3 which is the value obtained from QCD consideration [16] and ignore possible isospin breaking effects and the strange quark content in the nucleon. These can be incorporated as given in [17]. As can be seen above direct detection can strongly constrain λ ΦH

A. Scan strategy
A numerical study of the parameter space is performed as follow. We scan the full parameter space according to the following order: • A value of T dec is randomly chosen in the range between m µ and 2 GeV.
The lower bound is chosen to avoid the stringent experimental bound on K → π + (nothing). Furthermore, the phenomenology of scalars as light as that was discussed in [1] and we will not repeat it. Once T dec and M S are fixed, λ SH is determined by Eq. (12).
The upper bound θ max (M S ) is given by the OPAL direct search for M S > 1 GeV [14] and the indirect bound from LHC run-I [15] for M S > 60 GeV as discussed in previous section. We only found a few viable solution for M S < 2 GeV in our numerical scan, however we did not exclude this possibility in our study. We set the upper bound of |θ| < 2 × 10 −3 for M S < 2 GeV which comes mainly from the rare B decays [18] which is much more stringent than the OPAL bound.
In our numerical scan we found no solution for DM lighter than 0.5 TeV. This can be understood as follows. Since the requirement of RGE improvement of scalar stability will lead to large scalar couplings. Roughly speaking, the larger scalar couplings the larger DM annihilation cross section, and hence the smaller relic density. So one needs heavier DM to lower this cross section in order to get the relic density in the right ballpark. On the other hand, for the same couplings, the heavier ρ gives higher relic density at freeze out. Hence; there is an upper bound on M ρ so that relic density is not so high as to over close the universe. This is conservatively chosen to be 4 TeV.
• With the above set of parameters generated we calculate the following With these parameters we check Γ inv H to make sure the sum of all the invisible decay channels is still smaller than the experimental limit. If so this parameter set will be accepted as viable solutions.
The lower bound is from the positivity of the scalar potential and the upper bound is from the perturbativity.
Here, a consistency check is made so that κ =κ − λ φS v S < 0. This ensures ρ is the DM candidate.
With all the above parameters fixed, except λ φ which has no low energy constraint, we can go on to check whether the relic density and the DM direct search bound [19] are both met. Otherwise, the process will start over again.
• Lastly, we randomly scan λ φ ∈ [0, 4π]. Since there is no known constraint we use the RGE to determine its viable value. 23116. If a Landau pole is encountered, or λ φ,S become negative, or any of the 1 At 1-loop level the running of y t is the SM one. It is well known, see [9], that µ SM V S is very sensitive to the initial value of y t (or m t ) for RGE running. Since we use the SM as the reference point and all we require is that the lepton number violation scale is below µ V S . The top quark mass uncertainty does not enter to affect our study and conclusions.
positivity conditions is violated, i.e λ ij < −2 λ i λ j , in the stability region of λ H the parameter set is discarded. We denote the scale where the vacuum instability happens , the set is considered viable 2 . If after some large number (10 5 ) of tries without success, the whole set of parameters will be discarded and the scan goes to first step again; otherwise we register the parameter set as one viable configuration.
In fact, M N is also a free parameter in our model. However, our numerical experiment could not find any viable solution for M N < 0.5 TeV and the numerical results are not very sensitive to the actual value when M N ∼ few TeV. Therefore, in our study we just set M N = 1 TeV as a benchmark.

B. Overview of the numerical results
The viable configurations are easy to get if one only requires that the scale of vacuum instability is higher than the SM one, . For later use, we define to quantify how much the improvement of the vacuum stability scale comparing to the SM case. To emphasize how the vacuum stability and RGEs affect the parameters in our model, here we focus on those configurations with R V S > 2 and we have generated 4000 such viable sets of parameters. This choice is arbitrary and we intend it for illustration purpose only. We found many configurations with R V S > 2 and the largest R V S we got is ∼ 11 using the scan algorithm just stated. This is in agreement with the expectation that singlet scalars or Higgs portal models tend to improve electroweak vacuum stability. With this algorithm and the computing resource at hands, we did not find the configurations where the scale is pushed all the way to the Planck mass.
To demonstrate this we display the details of two typical configurations:  The other features of our numerical results can be summarized as follow: • It is easier to find solutions when T dec > ∼ 1.3 GeV and M S , V S , κ are not very sensitive to T dec , Fig.3(a-c). Moreover, R V S does not seem to depend on T dec . So we will focus instead the parameters dependance on M ρ .
• Solutions show that M ρ is in between roughly 1.5−4 TeV and center at around 2.5TeV with larger R V S , Fig.3(d-f). Although the range that M ρ ∈ {0.5, 4}TeV is scanned in our numerical study, we found no solutions with M ρ < ∼ 1.5 TeV.
• The upper bound of λ Φ depends on λ ΦS near λ ΦS > ∼ 1.0, Fig.4(g). less, λ 2 SH < 10 −3 even for M s ∼ 100 GeV, see Eq. (12). Therefore, λ ΦH plays the leading role of improving λ H stability. By linear extrapolation, the beta function needs an extra move up the Higgs scalar potential stability limit at µ SM V S to 100 × µ SM V S That amounts to λ ΦH ∼ ± √ 2.4 ∼ O(±1). However, the negative solution is not viable. Because the sizable negative λ ΦH will quickly drive the λ HS into large negative value during the RGE running such that λ HS < −2 √ λ H λ S and violates a vacuum stability condition. This estimate agrees with the feature we found in the numerical study that λ ΦH centers at around +0.5 and it is not very sensitive to other parameters.  This boundary λ ΦS > ∼ 4Y 2 S can be clearly seen in our numerical result, Fig.5(a). In addition to the issue of vacuum stability, the scalar sector in our model also introduces the Landau pole problem. As already mentioned, the beta function of λ Φ is always positive and it could leads to Landau pole. Assuming that λ ΦH and λ ΦS are more or less constant during the RGE running between t 0 and t, then Eq.(16b) admits an exact solution for λ Φ : where α 1 = 10/(16π 2 ), α 2 = (λ 2 ΦH + λ 2 ΦS /2)/(16π 2 ), and α 3 = tan −1 (λ Φ (t 0 ) α 1 /α 2 ). The Landau pole appears at t Landau when the argument inside tangent becomes π/2. Or, And it is clear now why this model prefers a heavy DM ( > ∼ 1.4 TeV ) after taking into account the RGE running and the issue of vacuum stability. When DM is relatively light, close to 1.4 TeV, the total cross section is saturated by the channels with SM final states. Since |λ ΦH | is not sensitive to M ρ , we immediately expect that σv W/Z/H / σv total is inversely proportional to DM mass squared, as it is shown in Fig.6

A. Extra light scalar S and its mixing with the Higgs boson
As mentioned in Section II, all signal strengths take a universal value of cos 2 θ in our model due to the H − S mixing. In Fig.7(a), the correlation between sin 2 θ and M S from our numerical study is displayed as well the expected sensitivity of sin 2 θ by improving the signal strength measurements at the LHC14 with 3ab −1 luminosity. The parameter space with M S > ∼ 40 GeV or equivalently the large mixing angle in our model will be covered by LHC14. If no detectable deviation is found, this part of parameter space will be discarded.
On the other hand, if this large mixing region is not excluded by LHC14, the same parameter space can be further directly probed by future facilities as we discuss next.
One immediate consequence of the existence of a neutral scalar of mass few tens of GeV and sizable mixing is that the triple SM Higgs coupling, λ SM HHH = 6λ H v H = 3M 2 H /v H , will be reduced. The tree level triple Higgs coupling is given in Eq.(18) as λ HHH . In Fig.7(b), the deviation δ HHH = (λ HHH − λ SM HHH )/λ SM HHH is displayed. The deviation can be as large as 20% when M S ∼ 100 GeV and center around few percents for configurations with better RGE improvement. This Higgs triple coupling is expected to be probed to 50% at LHC14 with 3ab −1 luminosity [20] and ∼ 10% at CEPC [21]. So some of the parameter space with M S > 40 GeV in this model can be probed at future colliders.
Similarly, the quartic coupling of the SM Higgs will be modified from λ SM could reach ∼ 30% when M S ∼ 100 GeV, see Fig.7(c). Note that for small λ S , δ HHH ∼ (c 3 θ − 1) and δ 4H ∼ (c 4 θ − 1), therefore |δ 4H | > |δ HHH |. In principle the quartic coupling could be probed through the triple Higgs production but this cross section is hopelessly small to be searched for at any foreseen future facility.

B. Decays of S
Since now that the mass of S is in the range of few tens to one hundred GeV, more decay channels are opened up than in the case that M S < 1 GeV which was discussed in [1]. The various decay widths can be easily derived from the well known ones in the SM. The widths of the dominant decay channels are given at below: where S → bb and S → ωω take up about 90% of the total decay width for M S < M W , M Z .
There are also S → gg and S → γγ decays induced at the 1-loop level: where we only keep the 1-loop top quark contribution for S → gg decay. For S → γγ, the W −loop contribution dominates over the top-loop contribution and the two have opposite sign which give rise to O(1) loop factor which we neglect in both cases. However, the resulting branching ratio is smaller than 10 −2 (10 −4 ) for S → gg(γγ) and can be ignored.
When M S > M W , M Z , the 3-body decays S → W W * → W ff ′ and S → ZZ * → Zff are opened up and start to play a role. The total decay widths are where We have summed over all light final states and treat them as massless particles and and we excluded the S → W tb, S → Zbb, and S → Ztt modes since they are either kinematically forbidden or suppressed. Above the mass thresholds, the branching ratio for S → W W * (ZZ * ) can reach ∼ 10 −5 (10 −6 ) when M S ∼ 100 GeV and can be completely ignored.
Note that all the decay widths can be fully determined by a given set of θ, v S and M S . In all, the branching ratios and the total decay width are displayed in Fig.8. One can see that in most of the parameter space, S → ωω is the dominate decay channel which has the invisible final states. Even for those configurations with sizable S → bb decay branching ratios still suffer from the small production cross section and make the study of S a challenging task.  [21]. These machines will allow measurements of the properties of the SM gauge bosons at unprecedented precision. Since the SM gauge bosons couple to S though the S − H mixing, the Z → ff S and W → ff ′ S decays are now opened. These decays branching ratios are given by [16] Br(Z → Sff ) for Z → ff S where r Z = M S /M Z , and a similar expression for Br(W →Sff ′ ) Br(W →ff ′ ) by substituting g 2 / cos 2 θ W with g 2 also r Z → r W = M S /M W . Clearly, the mixing sin 2 θ plays a pivot role to determine the size of branching ratios. In Fig.9(a), we display the V → Sf f branching ratio normalized by V → ff modulated with the mixing. And the outcome of our numerical experiments are shown in Fig.9(b,c). For M S < 60GeV, our model predicts a branching ratio around 10 −8 − 10 −6 times of the SM Br(V → ff ). It is interesting that there is a lower bound for these branching ratios and this is understood as the scalar potential stability requires a relatively large S − H mixing.
In order to make the best use of the 3-body decays we note that the dominant branching modes S is into bb or ωω. In turn they lead to the signatures Z → ff +bb and Z → f +f + / E where / E denotes missing energy. The particularly interesting ones are f = b, µ, e. The invariant mass squared distribution, M 2 f f , is a very useful quantity for suppressing the SM background. Defining y f = where The kinematic lower bound can be safely taken to be zero even for y b . This distribution peaks at y near the kinematic limit due to the propagator effect since the charged fermion pair comes from a Z * . The dominant SM background is due to Z → f * f or ν * ν follow by the f * (ν * ) → f (ν) + Z * /W * /γ * with the virtual gauge boson going into the appropriate final fermions. The y distribution peaks at smaller values as seen in Fig.10. The SM branching ratios for Z → bb + / E is 5.25 × 10 −8 and for Z → µμ + / E is 1.07 × 10 −8 . The latter is a very clean signal to utilize. Furthermore, the SM background from Z → Z * h * is 10 −4 times smaller than the above and can be ignored.  In passing we also note that the decay Z → ωωνν will contribute to the Z invisible decays but only at level < 10 −6 . This will be difficult even at the Z-factory mode of the FCC-ee. For M S < M W the decay channels W → S + W * will also be open. The virtual W * will then decay into a fermion pair. The signals will be similar to the Z decays discussed before.
However, we find that it will not add any additional information. It also suffers from lower event rates at the FCC-ee compared to the Z.

D. DM bound state
Our solutions indicate that the DM ρ has mass in the TeV range. Furthermore the parameter λ Φ is mach smaller than λ ΦH and λ ΦS . For DM with a mass of a few TeV or higher, all the masses of other states except χ can be ignored. More importantly the dark scalar ρ can interact with each other through exchanging the relatively light S and H in the t-channel and this force is attractive. The relevant interaction is given by and sinceκ ≫ λ ΦH v H the s mediation dominates. As shown in Fig.5 the B ρ propagator.
We start with the ρρB ρ coupling vertex. If one writes the effective coupling between the bound state B ρ and ρ as By dimensional analysis, α B can be estimated to be α B ∼ (κ 2 /M ρ ).
Next, the decay width of B ρ is proportional to its wave function absolute squared at the origin times the decay amplitude squared Γ B ∝ |ψ(0)| 2 ×|M Bρ | 2 . The probability density for two ρ's to meet is |ψ(0)| 2 ∼κ 6 /M 3 ρ by dimension analysis. We rescale the decay amplitude square to make it dimensionless and it can be further broken into where the subscripts label the decay final state. By setting all final states massless, and with the help of Eq.(17), we immediately have 3 where we have dropped terms suppressed by O(v H /M ρ ). Since both γ sH = O(v 2 H /M 2 ρ ) and γ ff = O(m 2 f /M 2 ρ ) thus can be neglected, and The dimensionless factor |M Bρ | 2 in our numerical analysis is ∼ O(1) and the bound state decay width can be estimated : Finally, we put Γ B into the propagator squared and the annihilation cross section due to the B ρ resonant can be estimated to be where the factor Γ B /M B is inserted to take care the nearly on-shell B ρ decay.
When v ≪ 1, s ∼ M 2 B and there is almost no temperature dependence, we have and is the boost factor for indirect DM detection.
For a typical value that |M Bρ | 2 ∼ O(10 0 ) andκ ∼ 0.1M ρ we have the boost factor around 100. In our numerical study, we found that the branching ratio of DM pair annihilate into mono-energetic Majoron pair is a few to 40%, Fig.6(c). The boost factor will make σv (DM + DM → ωω) ∼ 10 −26 − 10 −24 (cm 3 /s). The sizable annihilation cross section opens up a possibility that the Goldstone bosons could be a component of the 'apparent' neutrino flux at E ν = M ρ in IceCube and other neutrino observatories. Moreover, they give rise to shower events and no tracks and it is mostly originated from the Galactic center.
In additional to the mono-energetic Majoron line, there is sizable fraction that DM pair annihilate into SS pair, see Fig.6(b). The mono-energetic S then subsequently decays into bb and ωω. Our numerical experiment indicates that Br(S → ωω) > Br(S → bb) in most of the parameter space, Fig.8 41) is less important than the T ρ contribution. Therefore the σv becomes: Comparing to the σv 0 without BS, which is mainly controlled by Eq. (17), we gain an enhancement factor about for T ρ ∼ M ρ /20 andκ < ∼ 0.4M ρ . So that the bound state effect at the DM freeze out era is not important comparing to the direct DM-DM tree-level annihilation in our model.

E. Kinetic decoupling between DM and DR
Even after the thermal decoupling between the DM(ρ) and DR (ω), they can still interact with each other through the scalar quartic coupling term 1 4 λ ΦS ω 2 ρ 2 . Given that λ ΦS is sizable in this model, we would like to know whether this has any detectable cosmological implication. And the relevant question to ask is at what temperature, T k , the two will decouple kinetically.
By straightforward calculation, one has the nonrelativistic cross section for ω + ρ → ω + ρ for Majoron energy is much less than M ρ . After taking the thermal average, the rate for a single DM particle to collide with a Majoron is given by where n ω is the DR number density which behaves like that of photon and T is the temperature 4 . At low temperatures, the typical momentum of ω, p ω ∼ O(T ), is much less than the typical momentum of DM, p ρ ∼ O( M ρ T ). Therefore, for a DM to acquire a momentum transfer which is comparable to p ρ , it needs to accumulate many tiny momentum transfers from multiple collisions with the ambient DR. This process is very similar to the random walk and the number of collisions can be estimated to be √ N col p ω ∼ p ρ or N col ∼ M ρ /T .
And the kinetic decoupling temperature can be estimated by requiring that or Using the above estimate, we obtain the corresponding T k = 0.67(2.34) MeV for configuration- Thus, the temperature T k determines a lower bound on the masses of the smallest halos from the Jeans mass [23], Currently, the highest kinetic decoupling temperature can be probed is around 10keV, and the T k in our model is too high to be detected with the current observational precision.

V. CONCLUSION
We calculated the 1-loop beta functions for the minimal singlet Majoronic model[1] and performed a thorough numerical study on the parameter space of this model. In order to have an operational type-I see-saw mechanism, it is required that the lepton number breaking scale, v S , is lower than the scale µ V S where the SM electroweak vacuum become unstable.
The extra scalar degrees of freedom always help to improve the stability of SM electroweak vacuum, thus µ V S > µ SM V S . However, the right-handed Majorana neutrinos contribute negatively to the beta function for λ S through the Yukawa Y S . Additional attention to this new instability has to be taking into account and ensure that λ S > 0 when energy scale µ < µ V S . Moreover, the beta function for λ Φ is always positive so we looked for the solutions that there is no Landau pole below µ SM (1−loop) V S . Other phenomenological requirements had been considered in our numerical scan are: (1) the upper limit of SM Higgs invisible decay width, (2) Majoron decouple from the thermal bath at the temperature between m µ and 2 GeV, (3) the upper limit on the mixing between the SM Higgs and the beyond SM scalar, (4) the correct DM relic density, (5) the upper limit of direct DM searches.
The results of our numerical experiment have been summarized and discussed in Sec.III.
Here we highlight the physics of our finding. 1. A decoupling temperature at or below 2GeV leads to small λ SH , Eq.(12).
2. For the sake of λ S -vacuum stability, Eq.(16c), y S = √ 2M N /v S cannot be too large.
This leads to v S in the 2 − 20 TeV range which is relative large compared to the SM VEV.
3. In order to have such a value for v S , a large mixing angle between S and the SM Higgs is preferred, Eq.(22).