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 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.


Introduction
The strong first order phase transitions (PTs) are of great cosmological significance as they are necessary conditions for successful baryogenesis mechanism [1]-which explains why there are more matter than anti-matter in our universe. Since the first discovery of gravitational waves (GWs) by the Advanced Laser Interferometer Gravitational Wave Observatory (aLIGO) [2], a new era of GW astronomy has been opened, we are now at the position of exploring the multi-band GW detections. On the other hand, in literature it is well-studied that stochastic GW backgrounds can be generated during first order cosmological phase transitions [3][4][5][6]. Unfortunately, the electroweak phase transition (EWPT) in the standard model (SM) of particle physics turns to be a crossover [7][8][9] , which is far from strong enough to generate the stochastic GW backgrounds. In order to strengthen the EWPT as well as probe new physics, many extensions of the SM have been studied, such as models with extra scalar singlet(s) [10][11][12][13][14], with dimensional six effective operators [15][16][17][18][19], 2HDMs [20][21][22][23], NMSSM [24][25][26], and more other works [27][28][29][30].
The singlet Majoron model is another simplest extension of SM [31] (see also the GR model [32]), which is originally proposed to solve the neutrino mass problem. With the introductions of one more complex Higgs singlet conserves the global U (1) B−L symmetry and also the right-handed (RH) neutrinos, people get interested if there exist first order phase transitions [33][34][35][36] or testable GW signals in the developing GW astronomy as well as radio astronomy [37]. In Refs. [33,34] the authors claim there exists a flat direction on the surface of the effective potential, thus the two-field problem can be reduced into a single-field one, they conclude strong first order phase transitions can be realized. Another study in Ref. [35] further confirm that phase transitions 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 these earlier works, the researchers haven't considered the Yukawa couplings between the singlet Higgs field and the right-handed neutrinos. Later on, a very comprehensive work has been done in Ref. [36], the authors conduct a full numerical simulation to handle the two-field problem directly. Their work shows the existence of two step phase transitions-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 and the Yukawa coupling between the Higgs singlet and the right-handed neutrinos are expected in order to obtain strong first order phase transitions.
In this paper, we will follow the work of Ref. [36], treat the two-field problem directly by making use of the public available Python package CosmoTransition [38], and go even further by calculating the GWs generated from strong first order phase transitions. We confirm the patterns of phase transitions which have been found in Refs. [35,36], at the meanwhile we find even more new patterns, for the PTGWs we confront them with spaceborne gravitational wave interferometers such as LISA [39], DECIGO [40][41][42], BBO [43], and also the Chinese projects TAIJI (ALIA descoped) [44] and TianQin [45,46], we conclude they can be astronomically interested only when the vevs v BL have small values (smaller than ∼ O(10) GeV). The paper is organized as follows: in section 2 we write down the model of singlet Majoron and unify our notations; in section 3 we calculate the effective potential of the two classic fields within the finite temperature field theory, and show the interesting patterns of cosmological phase transitions; in section 4 we calculate the GW signals generated from first order phase transitions and confront them with space-borne GW detectors; finally in section 5 we make the conclusion.

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 apart from 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 doublets of SM Higgs and Light-handed(LH) fermions, f, 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 (v BL , v ew = 246 GeV), we will choose µ 2 To be clearer, let's expand the Higgs Bosons into their real components Just like a Dirac mass with m = f v ew /(2 √ 2) being obtained after the electroweak symmetry breaking, a Majorana mass with M = gv BL / √ 2 can also be generated after the breaking of the U (1) B−L symmetry. If the Majorana mass is much larger than the Dirac one, redefine heavy and light neutrinos as , then the masses of light neutrinos m ν m 2 /M are highly suppressed by the heavy neutrinos mass m N = M , that's the simplest way to understand why neutrinos' observational masses are so small.
A massless Goldstone boson χ appears after the U (1) B−L global symmetry broken, this is the Majoron field discussed in this paper. Even though the Majoron is massless, small mass term will be generated after one has considered the effective operators with higher dimensions. The authors of Ref. [47]  In [49,50], it is pointed a v BL with a value larger than v ew can be quite dangerous. However, in our numerical simulation we will explore both the scenarios either v BL be larger or smaller than v ew , aiming at getting a global overlook on the parameter space which can lead to strong first order PTs and GWs, remarkly we find v BL with higher values are also astronomically uninterested according to the results in Sec.4.

Patterns of cosmological phase transitions in the singlet Majoron model
Consider two classic fields the tree level potential can be rewritten as whose global minimum clearly shows today's vacuum. Then consider the one loop correction at zero temperature with MS renormalization with (c i , c W T , c Z T ) = (1.5, 0.5, 0.5), here 'T' stands for the transverse modes, 'L' stands for the longitudinal modes, the degree of freedoms (n H , n G , n ρ , n χ , n ν R , n t ) = (1, 3, 1, 1, −6, −12), (n W T , n W L , n Z T , n Z L ) = (4, 2, 2, 1), the mass spectrum can be found in the appendix. We also introduce a counter term ∆V T=0 ct (s, h) = Ah 2 , use the following renormalization conditions to make sure the vacuums at zero temperature won't be shifted away, we then solve and substitute the energy scale Q and A, see the appendix. The temperature correction reads here plus sign is for fermions and minus sign for bosons. In order to include the ring (or daisy) contribution, we make the following replacements [11,20] for 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 the appendix. We substitute the Debye mass terms into Eq.(3.6), finally, the total effective potential of the classic fields is given by Obtained the effective potential, bubbles' profile 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, it can be solved with a 'overshooting-undershooting' method [51,52], however the case for multi fields can be much more complicated, in the reference [38] Carroll L. Wainwright creates a powerful method-path deformation, to overcome the problem. In our project we make use of Wainwright's public available package CosmoTransition [38] to study the cosmological phase transitions. The strength of phase transitions can be illustrated by the parameter α, which is the ratio between the latent heat and the radiation energy where the latent heat is evaluated from ρ vac = [T (d∆V eff /dT ) − ∆V eff ] | T =T * , ∆V eff is the potential difference between the true and false vacuums; the radiation energy density is ρ r = g * π 2 T 4 * /30 with g * being the relativistic degree of freedoms. The time scale of phase transitions is the inverse of the parameter β where the Euclid action In the actual calculation, a renormalization of β could be quite usefulβ By definition, a larger value of α means a stronger phase transition, and a larger β means a faster phase transition. The phase transition temperature T * is estimated when the bubble nucleation probability in unit time and unit volume 1,  Fig. 1. An example of a first order phase transition along the h-direction followed by a higher order phase transition along the s-direction in the contour plot of the effective potential Eq. (3.8). This is given by the parameters v BL = 5 GeV, λ s = 1, λ sh = 0.1, g = 0.6. Here the h-direction phase transition happen at temperature T 137.88 GeV, the s-direction phase transition happen at the temperature T 10 GeV. this gives us a simple criterion to estimate T * (3.14) Notice when we derive this criterion, we have assumed the energy density of the universe is dominated by radiation, so it is not gonna be applicable when α is very large-which means the radiation dominated condition is spoiled-such as the case when bubbles have runaway in vacuum [53].
In the singlet Majoron model, parameters can be divided into two distinguished groups: the SM ones and the beyond-SM ones. For SM parameters, we take the gauge couplings g 1 = 0.65 and g 2 = 0.35, the heaviest fermion top quark's Yukawa coupling constant y t = 0.989, the Doublet Higgs boson's self coupling constant λ H = 0.13, today's electroweak vacuum expected value(vev) v ew = 246 GeV. The beyond-SM parameters are the newly introduced singlet Higgs and the right-handed neutrinos' self coupling λ s , mixing coupling λ sh , Dirac mass term Yukawa coupling f , Majorana mass term Yukawa coupling g, and vev of the global symmetry v BL . It easily to find the parameter f is degenerated with g and v BL , we thus won't set it as a free parameter. Finally, we have four parameters to be constrained: λ s , λ sh , g, v BL . For the plotting convenience, we will use a numeric string such as '100-0.8-0.1-0.4' to represent v BL = 100 GeV, λ s = 0.8, λ sh = 0.1, g = 0.4. In this paper we will scan the following parameter space (3.15) Since there are two scalar fields, the effective potential Eq.(3.8) is a 2D surface evolves with the changing of temperature. At early universe the temperature is high, symmetries are unbroken, the zero point of the fields space turns to be the only vacuum. As the universe expanding its temperature goes down, structure of the effective potential Eq.(3.8) becomes nontrivial and different local minimums will be formed, among them the lowest one is the true vacuum. Phase transitions can happen by tunneling from the higher vacuum(s) to lower ones, sometimes there are barriers between the distinguished vacuums which corresponds to the first order phase transitions, while at some other times phase transitions can happen very smoothly without any potential barriers being formed, they are referred as higher order phase transitions. In our numerical simulation, we find patterns of phase transitions 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, and the meanwhile first and/or higher PTs may occurs under the same parameters. Since we are mostly interested in the strong first order PTs, we will show some examples in the figure plots.
In Fig.1 we show an example of a first order phase transition along the h-direction followed by a higher order phase transition along the s-direction. As the universe cools near to the phase transition temperature T 137.88 GeV, there exist multi-vacuums and the tunneling happen along the Higgs direction firstly; another smoother phase transition happen when the temperature goes down to T 10 GeV, however it can not be seen in Fig.1 due the smallness of v BL . In this example the h-direction first order phase transition turns to be to too weak to generate detectable GW signals.
In Fig.2 a higher order PT happen around T 175 GeV along the s-direction, followed by a strong first order PT at the temperature T 86.6 GeV. We see as temperature goes down, the original vacuum located at the origin point splitting into two symmetric new vacuums on the axes of Higgs by a higher order phase transition. When the temperature goes even lower, two new vacuums on the h axes are formed, they coexist with the old vacuums until the tunneling happen at T 86.6 GeV. The first order PT is very strong, as a matter of  Fig. 4. An example of two first order phase transitions, the first is along the s-direction followed by the second which is along a mixing direction. This is given by parameters v BL = 100 GeV, λ s = 0.3, λ sh = 0.2, g = 0.3. Here the s-direction phase transition happen at temperature T 364.07 GeV, the hybrid phase transition along the mixing direction happen at T 122.7 GeV.
fact the GWs generated can be detected by the U-DECIGO detector [54], see the first plot in Fig.5. This strong first order phase transition may be regarded as the generalization analogy to the EWPT in the standard model. However, they are somehow different due to the mixing effect between Higgs and the s field, here we'd like to regard this case as an example of 'hybrid phase transition'. Fig.3 shows a first order phase transition along the s-direction followed by a higher order PT along the Higgs direction. This s-direction PT proceeds very fast, leads to very high frequency GWs which go beyond the detectable range of spatial GW interferometers. The final vacuums locate at (±v BL , ± v ew ) when the temperature gets closed to 0.
It is even more interesting to find multi first order PTs from some parameters, in Fig.4 we show an example with a h-direction first order PT happen after a s-direction first order PT. Near T 364.07 GeV, there form three vacuums separated by the barriers, one locating at the zero point, the other two are the nonzero ones with opposite signs according to the symmetry. As temperature goes lower, the true vacuum tunnels from the zero to the nonzero ones. Notice here we have employed zoomed-in plots to show the s-direction first PT better, as we have claimed, the PTGWs have too high frequency and too weak amplitude to be detectable. The other first order PT happen to be a hybrid one at T 122.7 GeV, whose GWs are also too weak to be detectable.
Therefore we see the cosmological phase transitions patterns in singlet Majoron model are quite fruitful, which should be the distinguishable feature for 2D problems compared with single field ones. We thus believe the single Majoron model is a useful laboratory to study the properties of cosmological PTs and PTGWs.

Phase transition gravitational waves in the singlet Majoron model
There are various GW sources generated during strong first order phase transitions, we will follow the paper Ref. [53] by considering the contributions from the scalar fields, the sound waves [55][56][57] , and the magnetohydrodynamic (MHD) turbulence [58]. The first GW source comes from the vacuum bubbles collision [6,58,59], the GW energy spectrum and peak frequency are this experience formulae are summarized from numerical simulations with envelope approximation (also see the work with analytic derivation [60,61] or beyond the envelope approximation [62,63]). The second source is generated by the sound waves of the bulk motion [55][56][57] Ω sw (f sw )h 2 2.65 × 10 −6 1 β The third are GWs generated by MHD turbulence [58] with energy spectrum and peak frequency The benchmark energy scales are taken v BL = 1 GeV, 5 GeV, 10 GeV, 50 GeV, 100 GeV, 246 GeV, with the notation 'v BL − λ s − λ sh − g' in the plots (here 'm0.3' refers to −0.3). The solid lines refer to GWs from h-direction or hybrid first order PTs, the dotted lines refer to those from s-direction PTs. For parameters which lead to twice first order PTs, we have use suffixes (h) or (s) to make it clearer.
where h * = 1.65 × 10 −5 Hz (T * /100 GeV) (g * /100) 1 6 . In Ref. [53] the authors discuss three different kinds of bubbles, i.e. the non-runaway bubbles, the runaway in plasma bubbles, and the runaway in vacuum bubbles. In this paper we will follow them directly. 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 when the bubble wall velocity is quite close to 1, or κ ν = v 6 5 w 6.9α (1.36 − 0.037 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 s field itself is also included, the smallest α in this case can be found by calculating where the squared mass difference between the two phases can be read from Eq.(A.1), the coefficients c a satisfy (c H , c G , c ρ , c χ , c W , c Z , c t , c ν R ) = (1, 3, 1, 1, 6, 3, 6, 3). The efficiency factor is evaluated from κ ν = (α ∞ /α)κ ∞ , with κ ∞ = α ∞ / 0.73 + 0.083 √ α ∞ + α ∞ .
When φ * /T * is very large, the phase transition may haven't happened until the universe is super-cooled, which could be quite different from the non-runaway or the runaway in plasma cases. In our paper, we will not deal with runaway in vacuum bubbles, since the criterion of T * Eq.(3.14) is spoiled when the total energy density of the universe is dominated by the vacuum energy. We scan the parameter space in Eq.(3.15) by splitting v BL into two groups: v BL ≤ v ew and v BL > v ew . As there are four free parameters v BL , λ s , λ sh , g, we will firstly fix the value of v BL as the benchmark energy scale, then increase the other parameters from their minimum to their maximum by a step length each time. For parameters which have the the same values of (v BL , λ s ), we pick out the ones which can produce the largest GW signals and display their GW signals in Fig.5 and Fig.6 by confronting them with the space-borne interferometers such as LISA, DECIGO, BBO, TAIJI and TianQin. Since there are so many curve lines in the plots we also take a unique notation 'v BL − λ s − λ sh − g' to save the plots space. Fig.5 shows the GWs from the singlet Majoron model with v BL ≤ v ew , where the benchmarks are chose to be v BL = 1 GeV, 5 GeV, 10 GeV, 50 GeV, 100 GeV, 246 GeV. It is clear that the detectable GW signals appears only when the v BL is small (smaller than ∼ O(10) GeV), on the other hand most of the GW curves are untouchable by interferometers since their amplitudes are not large enough while their peak frequencies are too high. By considering the different patterns of PTs discussed in Sec.3, we remarkably find the fact that all the detectable GWs belongs to the PTs along mixing directions, i.e. the hybrid phase transitions (HPTs). To get a better understanding about these HPTs, let's have a closer look at Fig.2, at the phase transition temperature T * 86.6 GeV the PT happens due to the tunneling from the false vacuum locates at 232 GeV, −5.43 × 10 −6 GeV to the true vacuum locates at 6.62 × 10 −6 GeV, 223 GeV , giving a φ * /T * √ 232 2 + 223 2 /86. 6 3.72 which undoubtedly shows a very strong first order phase transition. Other than HPTs, usually the first order PTs along h or s directions turns to produce GWs which are hardly detect by future interferometers. In our simulation, the strongest GW signals are given by the parameters v BL = 1, λ s = 0.1, λ sh = 0.2, g = 1.1, which has reached the detectivities of LISA, DECIGO, BBO and TAIJI. However, this sample may be quite closed to the parameters which can generate runaway in vacuum bubbles, as when we adjust parameters around it we continually meet very large value of φ * /T * (larger than ∼ O(10)) and spoilings of the criterion Eq.(3.14). Since we haven't dealt with the bubbles runaway in vacuum properly in this work, we may have lost some parameters which can generated detectable GW signals. Fig.6 shows GWs from vev v BL > v ew , where benchmarks energy scales are taken v BL = 500 GeV, 750 GeV, 1000 GeV. This time there exist none detectable signals, all the GWs are too faint while their frequencies are too high. And for larger a v BL , first order PTs along the s-direction tends to generate stronger GWs compared with the those from PTs along h-direction.

Conclusion
In this paper, we have studied cosmological phase transitions and gravitational waves in the singlet Majoron model, by treating it as a two-field problem fully without freezing any one of them. We confirm the phase transition patterns properties have been found in Refs. [35,36], i.e. usually the EWPT happens after the global U (1) symmetry breaking PT, further more we also find other patterns with different orders or along distinguished directions. Thanks to the 2D structure of the effective potential, PTs in the singlet Majoron model are much more fruitful compared with the ones in single field problems.
We further study the GWs generated during first order phase transitions, by separating the parameters with v BL which is lower and greater than electroweak v ew , we have scanned the parameter space, calculate and confront the GW signals with future space-borne interferometers such as LISA, DECIGO, BBO, TAIJI and TianQin. We find there could be astronomically interested GW signals, most of them are coming from first order phase transitions along the Higgs direction or a mixing direction, the space-borne GW interferometers may have chance to touch the parameters with small v BL such as smaller than ∼ O(10) GeV. For parameters with larger v BL > v ew the GWs from s-direction PTs become the domination source, however, they are too faint and with too high frequency to be detected.
We'd like to point out the fact that in this paper we just consider the bubbles which are non-runaway or runaway in plasma in Ref. [53], while the runaway in vauum bubbles (see [64]) are abandoned, for whom the criterion Eq.(3.14) is spoiled. By doing so we may lost some parameters which may lead to detectable GWs, we will leave this for future's work. no higher dimensional effective terms. We can diagonalize the mass matrix of neutral bosons, and thus redefine the ρ particle and the Higgs particle