Two-field cosmological phase transitions and gravitational waves in the singlet Majoron model

In the singlet Majoron model, we study cosmological phase transitions (PTs) and their resulting gravitational waves (GWs), in the two-field phase space, without freezing any of the field directions. We first calculate the effective potential, at one loop and at finite temperature, of the Standard Model Higgs doublet together with one extra Higgs singlet. We make use of the public available Python package ‘CosmoTransitions’ to simulate the two-dimensional (2D) cosmological PTs and evaluate the gravitational waves generated by first-order PTs. With the full 2D simulation, we are able not only to confirm the PTs’ properties previously discussed in the literature, but also we find new patterns, such as strong first-order PTs tunneling from a vacuum located on one axis to another vacuum located on the second axis. The two-field phase space analysis presents a richer panel of cosmological PT patterns compared to analysis with a single-field approximation. The PTGW amplitudes turn out to be out of the reach for the space-borne gravitational wave interferometers such as LISA, DECIGO, BBO, TAIJI and TianQin when constraints from colliders physics are taken into account.

The singlet Majoron model is a simple extension of the SM [42] (see also the famous GR model [43]), which was originally proposed to solve the neutrino mass problem. In this model, an additional complex Higgs singlet is introduced, which breaks the global U (1) B−L symmetry with its VEV and generates mass for the right-handed (RH) neutrino. Recently, there has arisen strong interest in the possibility of first-order phase transitions [44][45][46][47], GWs in radio astronomy [48] and dark matter physics [49]. In Refs. [44,45] the authors claim that a flat direction exists on the surface of the effective Majoron potential, allowing one to reduce the twofield problem to a one-effective-field problem. Their analysis shows that strong first-order PTs may be realized. Another study [46], further confirmed PTs can typically proceed in two steps, i.e., a very weakly first-order transition from the U (1) B−L breaking followed by the EWPT. In this earlier work, the Yukawa couplings between the singlet Higgs field and the right-handed neutrinos was not considered. Later on, a very comprehensive work was completed [47] by a full numerical simulation to handle the two-field problem. Their work shows the existence of two step PTs-and once again the PT due to the U (1) B−L symmetry breaking is above the EWPT. They also find large values of the two Higgs interactive coupling while the Yukawa coupling between the Higgs singlet and the right-handed neutrinos are required in order to obtain strong first-order PTs.
In this paper, we conduct a full numerical simulation to analyze the two-field cosmological PTs and evaluate the resulting GWs, by using the public available Python package CosmoTransition [50]. Thanks to the full 2D simulation, we are able to not only confirm the PT patterns which have been found in [46,47], but also find more new patterns, see Sect. 3 for more details. For each PT pattern we have found, we discuss the possibility to detect PTGWs with space-borne gravitational wave interferometers such as LISA [51], DECIGO [52][53][54][55], BBO [56], TAIJI (ALIA descoped) [57] and Tian-Qin [58,59].
The paper is organized as follows: in Sect. 2 we write down the singlet Majoron model and set up our notations; in Sect. 3 we calculate the effective potential of the two classic fields within the finite temperature field theory, and present patterns of cosmological PTs produced by our simulation that we classify; in Sect. 4 we evaluate the GWs generated from first-order PTs and put them in the light of future spaceborne GW detectors; finally, in Sect. 5 we summarize our main results and conclude.

The singlet Majoron model
The singlet Majoron model is among the simplest extensions of the Standard Model (SM) of particle physics. An U (1) B−L global symmetry, where the B and L stand for baryon and lepton numbers, is introduced in addition to the original gauge symmetry SU c (3) × SU (2) L × U (1) Y . Right-handed neutrinos ν R , together with a complex singlet Higgs σ become new members of the high energy physics particles zoo. The Lagrangian starts from in which Φ andL are the SM Higgs and Light-handed (LH) fermions doublets, and f and g are Yukawa coupling constants.
The potential for the doublet and singlet Higgs at tree level can be written as In order to ensure today's vacuum To be clearer, let us expand the Higgs Bosons into their real components Just as a Dirac mass with m = f v ew /(2 √ 2) is obtained after the electroweak symmetry breaking, a Majorana mass with M = gv B L / √ 2 is generated after the breaking of the U (1) B−L symmetry. If the Majorana mass is much larger than the Dirac one, we may redefine heavy and light neu- Then the masses of the light neutrinos, m ν m 2 /M, are highly suppressed by the heavy neutrinos mass m N = M. This is perhaps the simplest way to understand why neutrinos' observational masses are so small. However, this model cannot be the most satisfactory one when compared with the experimental data, from light neutrinos m ν = ( ; it would lead to an unnaturally small value for the Yukawa coupling f < 3 × 10 −7 √ g(v B L /GeV) 1/2 (by taking m ν < 2 eV). Therefore, in the present article we simply take it as a benchmark model in order to show in detail how multi-step phase transitions happen when there are more than one degrees of freedom being involved. The issue of unnatural smallness can be cured in some advanced models, namely, we refer to Refs. [60,61] for the associated analyses.
A massless Goldstone boson χ appears after the U (1) B−L global symmetry is broken, this is the Majoron field discussed in this paper. Although, with the above potential, the Majoron is massless, a small mass term can be generated after one has considered effective operators with higher dimensions. Following the arguments of Ref. [62], we adopt a constraint Since the neutrino mass upper bound is smaller than m ν ∼ 2 eV [63], it yields v B L 1.6 TeV.
In [64,65], it is pointed out that a v B L with a value larger than v ew can be dangerous. However, in our numerical simulation we will explore both scenarios with v B L larger or smaller than v ew , in order to obtain a global outlook on the parameter space which can lead to strong first-order PTs and GWs. Remarkably, we find v B L > v ew is also astronomically uninteresting according to the results in Sect. 4.

Patterns of cosmological phase transitions in the singlet Majoron model
Consider In vacuum at the present scale one gets back Eq. (3). At tree level the potential can be rewritten as whose global minimum clearly shows today's vacuum. Then we consider the one-loop correction at zero temperature with MS renormalization where (c i , c W , c Z ) = (3/2, 5/6, 5/6) and (n H , n G , n ρ , n χ , n W , n Z , n ν R , n t ) = (1, 3, 1, 1, 6, 3, −6, −12). The mass spectrum can be found in Appendix A. We also introduce a counter term ΔV T=0 ct (s, h) = Ah 2 and use the following renormalization conditions: to make sure the tree level vacuum at zero temperature is not shifted. We then solve and substitute the energy scale Q and A; see Appendix A. The temperature correction reads the plus sign is for fermions and the minus sign for bosons. In order to include the ring (or daisy) contribution, we make the following replacements [11,22,66] for the bosonic masses: The expressions for the self-energies and the thermal mass (Debye mass) for the longitudinal component of Z boson can be found in Appendix A. We substitute the Debye mass terms into Eq. (13). By doing so, the replacement Eq. (14) is applied in the effective potential (see the analyses in [22,66]), not only for the cubic term (see the analyses in [11,22,67]). The two different prescriptions can be understood as the theoretical uncertainty in one-loop perturbation theory, which can be safely handled. Finally, the total effective potential of classic fields is given by With the above effective potential, the bubble profiles can be found by solving the bounce equation where r 2 = |x| 2 , α = 2 at finite temperature, and r 2 = |x| 2 − t 2 , α = 3 at zero temperature. For single-field problems, they can be solved by the so-called 'overshootingundershooting' method [68,69]; however, the case of multiple fields becomes much more complicated. To overcome this situation, a powerful method dubbed 'path deformation' was developed in Ref. [50]. In our analysis, we make use of the public available package CosmoTransition [50] to study the PTs in the case of multiple fields. The strength of the PTs can be illustrated by the parameter α, which is the ratio between the latent heat and the radiation energy, where the latent heat is evaluated using ΔV eff is the potential difference between the true and false vacua; the radiation energy density is ρ r = g * π 2 T 4 * /30 with g * being the relativistic degrees of freedom. The time scale of the PT is the inverse of the parameter β, where the Euclid action is , and the bubble nucleation rate defined by Γ = Γ 0 exp (−S 3 /T ). In the actual calculation, a renormalization of β is quite useful: By definition, a larger value of α means a stronger PT, and a larger β means a faster PT. The PT temperature T * is estimated when the bubble nucleation probability in units of time and units of volume 1, This gives us a simple criterion to estimate T * : Note that on when we deriving this criterion, we have assumed that the energy density of the universe is dominated by radiation, so it is not applicable when α is very largewhich means the radiation dominated condition is spoiled, e.g., when bubbles show runaway in the vacuum [70].
In the singlet Majoron model, parameters can be divided into two distinct groups: the SM ones and the beyond-SM ones. For the SM parameters, we take the gauge couplings g 1 = 0.65 and g 2 = 0.35, the top quark's Yukawa coupling constant y t = 0.989, today's electroweak vacuum expected value (VEV) v ew = 246 GeV; the only underdetermined parameter is the doublet Higgs boson's self-coupling constant λ h , which may be shifted away from m 2 H,SM /(2 × v 2 ew ) due to the mixing effect. The beyond-SM parameters are the newly introduced singlet Higgs self-coupling λ s , the mixing coupling λ sh , the Dirac mass term Yukawa coupling f , the Majorana mass term Yukawa coupling g, and the VEV of the global symmetry v B L . It easy to see that the parameter f is degenerate with g and v B L once a light neutrino mass scale is imposed; we thus will not set it as a free parameter. Finally, we have five parameters to be constrained: v B L , λ s , λ h , λ sh , g. For convenience in plotting, we will use a numeric string such as '0.3-0.124-0.24-1' to represent For direct comparison of theoretical and experimental constraints in Refs. [71][72][73] according to the formulas in Appendix A-notice only three of them are independent after m 2 1 or m 2 2 takes the SM Higgs boson's mass. As a matter of fact, we find the most convenient path is as follows: firstly give values to v B L , λ s and λ h , then use their relationship to evaluate λ sh . When the SM Higgs is the lighter eigenstate, i.e., m 1 = 125 GeV, from Eq. (42) we have and in the opposite case, m 2 = 125 GeV, we can use to evaluate m 1 . In both cases In this research we scan the following parameter space: and by doing so we directly compare the parameters with the results in Ref. [73]. In Ref. [ 25 GeV v B L . In order to cover all the parameter space, we will employ the following scan strategy.
-Step I: start from the minima of v B L , λ s , and λ h . -Step II: evaluate the mass of the (non)SM Higgs from Eq. (22) or Eq. (23), λ sh from Eq. (24), sin θ and tan β from Appendix A. -Step III: compare (tan β, m 2 1 , m 2 2 , sin θ) with the theoretical and experimental constraints. If they can pass, then go to Step IV, otherwise iterate v B L , λ s and λ h and go to Step II.
starts form its initial value) into the model and calculate the PTs, iterate g.
Step II, or end when the all the parameter points have been exhausted.
Since there are two scalar fields, the effective potential of Eq. (15) is a 2D surface and evolves with changing temperature. In the early universe the temperature is high, symmetries are unbroken, the zero point of the fields' space turns out to be the only vacuum. As the universe expands its temperature falls, the structure of the effective potential of Eq. (15) becomes nontrivial and different local minima will be formed, among them the lowest one is the true vacuum. PTs can happen by tunneling from the higher vacuum or vacua to lower ones. Sometimes there are barriers between the distinct vacua, which correspond to first-order PTs, but at some other times PTs can happen very smoothly without any potential barriers being formed; these are referred as higher-order PTs. In our numerical simulation, we find that patterns of PTs in the singlet Majoron model can be quite fruitful, usually there are multi-step PTs, some of them are along the s field direction, some are along the h field direction, and some others are along a mixing direction. Since we are mostly interested in the strong first-order PTs, we will show some examples in the figures.
In Fig. 1 a higher-order PT occurs around T 182 GeV, along the s field direction, followed by a strong first-order PT at the temperature T 118.89 GeV. We see that as temperature falls, the original vacuum located at the origin point splits into two symmetric new vacua on the Higgs axis; this is a higher-order PT. When the temperature drops even further, two new vacua on the h axis are formed. These coexist with the old vacua until the tunneling occurs at T 118. 89 GeV. This first-order PT is very strong; here it tunnels from one vacuum located on the first axis directly to the second one located on another axis. Due to the mixing effect between Higgs and the s field, hereafter we shall refer Fig. 1 An example of a higher-order s-direction PT followed by a hybrid first-order PT. This is given by the parameters v B L = 25 GeV, λ s = 0.35, λ h = 0.126, λ sh = 0.388, g = 1.5. Here the higher-order PT happens around T 182 GeV, the hybrid PT happens at T 118.89 GeV to these as 'hybrid PTs'. In our simulation, we do find some parameters which can lead to very strong hybrid PTs and detectable GWs for the mentioned interferometers. Unfortunately, they also lead to unacceptable branch ratio for the H → hh decay channel. For the parameter in Fig. 1, it predicts m 1 = 7.9 GeV, and BR H →hh 0.99, which definitely has been ruled out by collider experiments. We present it here in order to completely show the 2D PT properties of the singlet Majoron model. Figure 2 shows a first-order PT along the s-direction (notice the zoomed-in plots) followed by a higher-order PT along the mixing direction. The occurrence of the s-direction PT can also be attributed to the mixing coupling λ sh , from Fig. 2 one can see s/T ∼ 100/357 < 1, which means the PT is weak. Indeed, it proceeds very fast, leads to very high frequency GWs which go beyond the detectable range of spatial GW interferometers. As mentioned before, the final vacuum is located at (v BL , v ew ) obviously when the temperature gets close to 0. This parameter gives m 2 328 GeV and thus belongs to the high mass region (see the details in Table 3).
It is even more interesting to find multiple first-order PTs from some parameters. In Fig. 3 we show an example with a hybrid first-order PT happening after a s-direction firstorder PT. Near T 294.06 GeV, there exist three vacua separated by barriers, one located at the zero point; the other two are the nonzero ones with opposite signs according to the symmetry. As temperature drops, the field tunnels from the origin to the nonzero VEV. Notice here we also employ zoomed-in plots to better show the two first PTs. As above the PTGWs have too high frequency and too weak amplitude to be detectable. The other first-order PT happens to be a hybrid one at T 132.23 GeV; its GWs are also too weak to be  Table 2).
We also find that, for some parameter space, it would allow weak first-order PTs above the U (1) symmetry breaking, but the perturbative analysis has become jeopardized if the Higgs VEV in the broken phase is too small. To address this issue, one probably needs to do a lattice analysis (namely following [9]) that would find a crossover in this parameter space. However, as we have seen from the above analyses, the patterns of cosmological PTs in a singlet Majoron model are already very fruitful, which is a distinct feature for 2D problems compared to single-field ones. We thus believe that the single Majoron model is a useful laboratory to study the properties of 2D cosmological PTs and GWs. In the followup study we plan to develop our numerical method further as regards the lattice analysis on the issue of non-perturbative regime in order to explore more details on cosmological PTs and GWs in higher dimensions.

Phase transition gravitational waves in the singlet Majoron model
There are various GW sources generated during strong firstorder PTs; we follow the study in [70] by considering the contributions from scalar fields, sound waves [74][75][76], and magnetohydrodynamic (MHD) turbulence [77]. The first GW source comes from the vacuum bubble collision [6,77,78], the GW energy spectrum and peak frequency are Here, the s-direction PT happens at temperature T 294.06 GeV, the hybrid PT along the mixing direction happens at T 132. 23 GeV f env 1.65 × 10 −5 Hz 0.62 These formulas found by experience are summarized from numerical simulations with envelope approximation (also see the work with an analytic derivation [79,80] or beyond the envelope approximation [81,82]). The second source is generated by the sound waves of the bulk motion [74][75][76], f sw 1.9 × 10 −5 Hz 1 v w β T * 100 GeV g * 100 1 6 . (29) The third are GWs generated by MHD turbulence [77] with energy spectrum and peak frequency where h * = 1.65 × 10 −5 Hz (T * /100 GeV) (g * /100) 1 6 . In Ref. [70] the authors discuss three different kinds of bubbles, i.e., the non-runaway bubbles, the bubbles with runaway in plasma, and the bubbles with runaway in vacuum. In this paper we will follow them directly (also see Ref. [83] for more details). For non-runaway bubbles, we evaluate the total GW energy spectrum from Ω sw ( f sw )h 2 + Ω tu ( f tu )h 2 , take the efficiency factor κ ν = α(0.73 + 0.083 √ α + α) −1 when the bubble wall velocity is quite close to 1, or κ ν = v 6 5 w 6.9α(1.36 − 0.037 √ α + α) −1 when v w is smaller than 0.1; we take the turbulence efficiency factor κ tu = 0.05κ ν . For the case of runaway bubbles in plasma, the contributions from the s field itself is also included, and the smallest α in this case can be found by calculating where the squared mass difference between the two phases can be read in Appendix A. The coefficients c a satisfy 1, 3, 1, 1, 6, 3, 6, 3). The efficiency factor is evaluated from very large, the PT may not have happened until the universe is super-cooled, which is quite different from the non-runaway case or the case of runaway in plasma. In our paper, we will not deal with runaway in vacuum bubbles, since the criterion of T * of Eq. (21) is spoiled when the total energy density of the universe is dominated by the vacuum energy.
Carrying out our scan strategy, for parameters with the same (v B L , λ s ), we pick out the ones which can produce the largest GW amplitudes and display their signals in Fig. 4 with v B L < 246 GeV and Fig. 5 with v B L ≥ 246 GeV, meanwhile confronting them with the space-borne interferometers such as LISA, DECIGO, BBO, TAIJI and TianQin. Note that there are too many curves in some subplots; in order to save the plots space we will take the unique string notation such as 'λ s − λ h − λ sh − g'. For the benchmark parameters taken in Figs. 4 and 5, we also make three tables (Tables 1, 2 for v B L < 246 GeV and Table 3 for v B L ≥ 246 GeV) in Appendix B to show more details, especially as regards which mass region (HM, IM and LM) they belong to. We can use the three tables to analyze how they could survive from the theoretical and experimental constraints such as in Ref. [73].
When v B L < 100 GeV, the most of the benchmark parameters lead to (non-)SM Higgs with mass m 1 < m H,SM /2, and they thus are severely constrained by the collider experiments on the decay channel H → hh. Since the upper bound about the decay branch ratio BR H →hh 0.2 [84,85], the coupling λ h must be quite close to the SM value m 2 H,SM /(2 × v 2 ew ), and sin θ nearly equals 1. In order to survive from the H → hh constraints, in Table 1 Fig. 4, we show the largest GWs generated for the benchmark scales. All the display curves are due to the hybrid first-order PTs. Unluckily these GWs are hardly within reach of the planned interferometers.
For 100 GeV ≤ v B L < 246 GeV, all the three cases (LM, IM and HM) can be realized. In Table 2 we show the parameters which have been taken in Fig. 4 for benchmark scales v B L = 100 GeV, 150 GeV, 200 GeV. Taking advantage of Fig. 7 in Ref. [73], one can see that a large value (close to unity) for sin θ is needed for the LM cases, while a small value (smaller than ∼ 0.4) is needed for the HM cases. We have compared the parameters in Table 2 to make sure that they are not ruled out by the experiments. Notice again that there exists an upper bound for λ s , i.e., 0.8 for 100 GeV, 0.5 for 150 GeV, and 0.45 for 200 GeV. This is due to the constraint contour on tan β and m 2 in the HM region; see Fig. 10(b) in Ref. [73]. Compared with the smaller v B L in Table 1, there tend to be more and more GWs generated by multiple first-order PTs.
In Fig. 5 we show the GWs generated by parameters with v B L ≥ 246 GeV; see Table 3 for more details as regards the chosen parameters. All except one belong to the HM region, the upper bound for λ s becomes smaller and drops slower. As we have mentioned, for these high v B L parameters, the PTs along the s-direction play a more and more important role in generating GWs, with even higher peak frequencies.

Conclusion
In this paper, we study the properties of cosmological PTs, especially the resulting GWs, in the singlet Majoron model, using a numerical simulation to treat the two-field problem without freezing any of the field directions. Compared with a single-field treatment, the patterns of PTs turn out to be much more diverse. We have not only verified the pattern with an EWPT happening after the global U (1) symmetry breaking, but also we find new patterns, such as strong hybrid PTs happening before the U (1) symmetry breaking. Our simulation suggests that the single Majoron model is an ideal benchmark in understanding the phenomenology of two-field cosmological PTs.
The PTGWs are likely not within the reach of detectability of space-borne interferometers such as LISA, DECIGO, BBO, TAIJI and TianQin-either their amplitudes are too low or their frequencies are too high for those next-generation instruments. As a matter of fact, without considering the collider experimental exclusion bounds, we are able to find strong hybrid PTs which can generate detectable GWs. Unfortunately, in such cases, the mixing coupling necessarily has to be large, and the non-SM Higgs need acquire a mass smaller than the half mass of the heavier Higgs. Experimentally this parameter space is ruled out by constraints on the H → hh decay branch ratio. However, we emphasize that the aforementioned conclusion is model dependent and the numerical method developed in the present study shall be widely applied to another model of cosmological PTs, which involves more than one degrees of freedom.
Finally, we would like to highlight the implications of the developed numerical method, which could be extended from several perspectives in future study. Note that we have only considered the bubbles which belong to the types of non-runaway or runaway in plasma (see Ref. [70]), but we abandon those bubbles with runaway in vacuum (see Ref. [86]; see also Ref. [87] for a most recent study), since they spoil our criterion in Eq. (21). By doing so we expect that the parameter points would be further constrained. Additionally, the same analysis may be applied to some cosmological PTs that occur at extremely low frequency band and hence might be falsifiable in the forthcoming experiments of primordial gravitational waves [88][89][90][91].
respectively. Today we have (s, h) = (v B L , v ew ), and hence, Here the subscript '(0)' refers to the results in the present vacuum. We can discard it without bringing about any ambiguity, and then today's mixing angle satisfies In order to compare with other work, let us also introduce another angle, Consider a counter term from the renormalization conditions of Eq. (11), we can get the energy scale cut-off and coefficient A, where (c i , c W , c Z ) = (3/2, 5/6, 5/6).
The self-energies are given by

Appendix B: The benchmark parameters
In the second part of the appendix, we show the benchmark parameters in Figs. 4 and 5, which are presented in Tables 1,  2, and 3. Apart from the model parameters, we also provide the results for the parameterization of tan β, m 2 1 , m 2 2 , sin θ , to directly compare them with theoretical and experimental constraints such as in Ref. [73]. Following [73], we have classified the parameter samples into three groups, which are the LM, IM, and HM regions, respectively.