Gravitational Wave Signatures from Domain Wall and Strong First-Order Phase Transitions in a Two Complex Scalar extension of the Standard Model

We consider a simple extension of Standard Model by adding two complex singlet scalars with a $\rm{U}\left(1\right)$ symmetry. A discrete $\mathcal{Z}_2 \times \mathcal{Z}^{\prime}_2$ symmetry is imposed in the model and the added scalars acquire a non zero vacuum expectation value (VEV) when the imposed symmetry is broken spontaneously. The real (CP even) parts of the complex scalars mix with the SM Higgs and give three physical mass eigenstates. One of these physical mass eigenstates is attributed to the SM like Higgs boson with mass 125.09 GeV. In the present scenario, domain walls are formed in the early Universe due to the breaking of discrete $\mathcal{Z}_2 \times \mathcal{Z}^{\prime}_2$ symmetry. In order to ensure the unstability of the domain wall this discrete symmetry is also explicitly broken by adding a bias potential to the Lagrangian. The unstable annihilating domain walls produce a significant amount of gravitational waves (GWs). In addition, we also explore the possibility of the production of GW emission from the strong first-order phase transition. We calculate the intensities and frequencies of each of such gravitational waves originating from two different phenomena of the early Universe namely annihilating domain walls and strong first-order phase transition. Finally, we investigate the observational signatures from these GWs at the future GW detectors such as ALIA, BBO, DECIGO, LISA, TianQin, Taiji, aLIGO, aLIGO+ and pulsar timing arrays such as SKA, IPTA, EPTA, PPTA, NANOGrav11 and NANOGrav12.5.


Introduction
The progress of the detection of gravitational wave (GW) event from a binary black hole merger confirmed by the Laser Interferometer Gravitational-Wave Observatory (LIGO) group [1,2] results in a remarkable enhancement in the physics of cosmology and astrophysics. The GWs from the early Universe can provide significant information related to the high energy physics phenomena in early Universe as the GWs after the production propagate without suffering any interaction and conserve almost all necessary physics information [3]. The production mechanisms of primordial GWs are associated with various cosmological sources such as inflationary quantum fluctuations [4], preheating after inflation [5], strong first-order phase transitions in the early Universe [6,7], topological defects of the domain walls, cosmic strings [8]- [10] etc. In this work, we explore two different production mechanisms of GW namely annihilation of cosmic domain walls (DWs) and strong first-order phase transition (SFOPT), in a simple two complex scalar extended Standard Model (SM).
Domain walls, the two-dimensional surface-like topological defects are originated when a discrete symmetry is spontaneously broken [11]. But generally the formation of domain walls appears to contradict some basic wisdom in cosmology [12] because energy density of the domain walls affects the total energy density of the Universe and leads to a rapid expansion of the Universe which is disfavoured by the present observational results [13]. One can solve this domain wall problem by considering the unstable domain walls which collapse before they would overclose the Universe [14,15]. In a theory, the unstability of domain walls can be established by considering that the discrete symmetry is approximate and it is broken explicitly [13]. After the formation of unstable domain walls they collide or annihilate and emit gravitational waves. Such GWs can contribute as a stochastic background of GW in the present Universe. Different types of particle physics models are considered in Refs. [13,16,17] for the production of GWs from domain walls.
On the other hand, GWs can also be emitted from strong first-order phase transitions. Initially, at high temperatures, the Universe is in a false vacuum state but gradually as the temperature of the Universe decreases as it expands and evolves, it shifts to the true vacuum state through the process of tunnelling. The process of tunnelling leads to the first-order phase transition via bubble nucleation. These bubbles could be of different sizes. The bubbles that are large enough for avoiding collapse will expand, collide and coalesce which distorting the spherical symmetry of the bubbles and eventually leads to phase transition and emission of GWs as a result. In this scenario, the GWs are produced via three mechanisms such as (i) bubble collisions [18]- [23], (ii) sound waves induced by the bubbles running through the cosmic plasma [24]- [27] and (iii) turbulence induced by the bubble expansions in the cosmic plasma [28]- [32].
Although the first-order electroweak (EW) phase transition (PT) would have been possible in the framework of SM of particle physics through the Higgs mechanism but with the observed Higgs mass of 125.09 GeV [34]- [37] the transition is a smooth cross-over and not a first-order one. However, in literature there are ample references where the authors have shown that SFOPT can be realised by simple extension of SM [38] - [57]. In the present work however, we primarily explore GW emission from the annihilation of domain walls by proposing a simple extension of SM with two complex scalars. In addition we also furnish the GW emission from SFOPT within the same framework of our proposed model. Domain wall is formed when a discrete symmetry is spontaneously broken. It appears at the boundary of two domains that is produced following such spontaneous breaking of discrete symmetry. In this regard, as mentioned, we extend the SM by adding two extra complex scalars and impose a discrete Z 2 × Z 2 symmetry to the model. The two complex scalars acquire non zero vacuum expectation values (VEVs) when the imposed discrete symmetry is spontaneously broken. The real parts of the complex scalars mix with the SM Higgs and give three physical mass eigenstates. One of the components of the physical scalars attributes to the SM like Higgs boson with mass 125.09 GeV. As the discrete symmetry in this framework is broken spontaneously, domain walls are formed which is made unstable by adopting a bias term that lifts the degenerate vacua. In Ref. [16] the authors have considered single scalar extension of SM with a Z 3 symmetry and subsequent production of GWs from domain wall. In this work however we approach the formation and collapse of domain walls by proposing two scalars (complex) extension of SM with a Z 2 × Z 2 symmetry. In addition we also demonstrate that by the extension of the two new complex scalars the effective degrees of freedom of the model increase by four units hence SFOPT can be possible in the model. In Refs. [39]- [57] the authors have proposed different particle physics models for the production of GWs from first-order phase transition but variance from their approaches in this work, we consider an extension of SM with two complex scalars in such a way that it can provide GW production from SFOPT as well as the GW production from unstable domain walls. We constrain the model parameters by using some theoretical and experimental constraints such as vacuum stability, perturbativity and the collider bounds. We choose five benchmark points (BPs) from the allowed model parameter space to calculate the GW intensity and peak frequency of GW for both the cases namely domain wall and SFOPT. We also investigate the detectability of such GWs at the future space-based detectors such as ALIA [58], BBO [59], DECIGO [60], LISA [61], TianQin [62], Taiji [63], ground-based detectors such as aLIGO [64], aLIGO+ [64] as well as low-frequency pulsar timing arrays (PTAs) [65] such as SKA [66]- [68], IPTA [69]- [72], EPTA [73]- [75], PPTA [76,77], NANOGrav11 [78]- [81] and NANOGrav12.5 [82,83].
The paper is organised as follows. In section 2, we briefly describe our two complex scalars extended SM and also derive the necessary relations between model parameters. In section 3, we discuss some bounds which we have used to constrain the model parameter space. The formation of unstable domain walls and the calculations of intensities of the resulting GW emission for this extended model are presented in section 4. In section 5, we furnish the finite temperature effective potential to study the first-order phase transition properties in the framework of our model. The possible production mechanisms of GWs from the strong first-order phase transitions are also explored in this section. Then in Section 6, we calculate the intensity and frequency of gravitational waves produced from both the domain walls and strong first-order phase transitions. In this section, we also discuss the detectability of such GWs at the pulsar timing arrays and future GW detectors. Finally, we summarise our work and give some concluding remarks in section 7.

The Model
In this work, we extend the Standard Model of particle physics by adding two complex scalar singlets S 1 and S 2 . A discrete Z 2 × Z 2 symmetry is imposed in the model where under Z 2 symmetry S 1 transforms as −S 1 and under Z 2 symmetry S 2 transforms as −S 2 . Thus under Z 2 × Z 2 symmetry S 1 and S 2 are (−1, 1) and (1, −1) respectively. The complex scalars S 1 and S 2 acquire non zero vacuum expectation values when the imposed Z 2 and Z 2 symmetries are broken spontaneously. The real parts (CP even) of the complex scalars mix with the SM Higgs and give three physical mass eigenstates. Here, we consider one of the physical mass eigenstates behaves as the SM Higgs boson with mass 125.09 GeV [84].
The most general renormalisable scalar potential under the Z 2 × Z 2 symmetry with U (1) symmetry can be written as where H, S 1 and S 2 are the SM Higgs doublet and the two complex scalar singlets respectively. In the potential (Eq. 1) we consider the mass parameters and all the couplings are real. After spontaneous breaking of SU (2) L × U (1) Y symmetry SM Higgs field H acquires a non zero VEV v =246.22 GeV and due to the spontaneous breaking of Z 2 and Z 2 symmetries the complex scalars S 1 and S 2 gain non zero VEVs v 1 and v 2 respectively (v 1 and v 2 are VEVs for S 1 and S 2 respectively). Therefore SM Higgs and the additional complex scalar fields can be represented as where h , s 1 , s 2 , χ 1 and χ 2 are the unphysical scalars. In Eq. 2, G + , G 0 are the charged and the neutral Goldstone bosons after the spontaneous symmetry breaking (SSB). We minimise the scalar potential (Eq. 1) using the following conditions (3) The minimisation conditions lead us to the following tadpole relations The real parts of the complex scalars in Eq. 2 mix with the neutral component of the SM Higgs doublet h. Considering the 3 × 3 mass matrix in the h − s 1 − s 2 basis we diagonalise the mass matrix by a unitary transformation U (θ 12 , θ 13 , θ 23 ) to obtain three physical mass eigenstates h 1 , h 2 and h 3 in terms of the old basis h − s 1 − s 2 and the three mixing angles θ 12 , θ 13 , θ 23 . The unitary transformation U (θ 12 , θ 13 , θ 23 ) has the following form where U (θ 12 , θ 13 , θ 23 ) is the standard Pontecorvo-Maki-Nakagawa-Sakata (PMNS) [85] matrix with the three mixing angles θ 12 , θ 13 , θ 23 and complex phase δ = 0. The matrix U can be written as cos θ 12 cos θ 13 sin θ 12 cos θ 13 sin θ 13 − sin θ 12 cos θ 23 − sin θ 13 sin θ 23 cos θ 12 cos θ 23 cos θ 12 − sin θ 12 sin θ 13 sin θ 23 sin θ 23 cos θ 13 sin θ 12 sin θ 23 − sin θ 13 cos θ 12 cos θ 23 − sin θ 23 cos θ 12 − sin θ 12 sin θ 13 cos θ 23 cos θ 13 cos θ 23 so that Eq. 16 takes the form The three physical mass eigenstates h 1 , h 2 and h 3 are therefore In this work we consider h 1 to be the SM like Higgs boson with mass m h 1 = 125.09 GeV. Using the minimisation conditions in Eq. 4, one can obtain the following relations The model parameters (µ H , µ S 1 , µ S 2 , λ H , λ S 1 , λ S 2 , λ HS 1 , λ HS 2 , λ S 1 S 2 ) in Eq. 1 can be expressed in terms of the physical masses of the scalar particles m h 1 , m h 2 , m h 3 , VEVs of the scalar particles v, v 1 , v 2 and by the three mixing angles θ 12 , θ 13 and θ 23 . Thus we take m h 1 , m h 2 , m h 3 , v, v 1 , v 2 , θ 12 , θ 13 , θ 23 , g S 1 S 2 , and κ S 1 S 2 as the input parameters of the model and the benchmark points (BPs) are chosen for different sets of values of these parameters that satisfy all the constraints both theoretical and experimental described in the next section to explore the production of the gravitational waves from domain walls and strong first-order phase transition.

Constraints
In this section, we discuss some theoretical and experimental bounds which we have used to constrain the model parameter space.

Collider Constraints
In the present scenario, our model is extended by two new extra complex scalars and we expect that the newly added particles can affect the LHC collider physics phenomenology. In this work, we consider h 1 as the SM like Higgs boson with mass 125.09 GeV so we use the collider bounds to further constrain the model parameters. The signal strength of the SM like Higgs h 1 is given by [85] where Γ SM is the total SM Higgs decay width. In Eq. 25, Γ 1 denotes the total decay width of SM like Higgs boson of mass 125.09 GeV which has the following form where Γ inv 1 refers to the invisible Higgs decay width. Here two invisible decay channels for h 1 can be possible. One of them is Γ inv The total invisible Higgs decay width can be expressed by the sum of the decay widths of these two decay channels as The expressions for the decay channels are and where the couplings g h 1 h 2 h 2 and g h 1 h 3 h 3 can be expressed in terms of the mixing angles, VEVs and the model parameters. We compute these two couplings numerically. The invisible decay branching fraction for SM like scalar h 1 can be written as We adopt the bounds on the signal strength of SM Higgs R 1 ≥ 0.84 [87,88] and on the invisible decay branching fraction for SM Higgs, Br inv ≤ 0.24 [89] to further constrain the model parameters.

Formation of Domain Walls and the consequent production of Gravitational Waves
As mentioned earlier, a discrete symmetry Z 2 × Z 2 is imposed on the potential in (Eq. 1). When Z 2 and Z 2 are spontaneously broken the fields S 1 and S 2 acquire non zero VEVs and domain walls are formed around the boundaries of the newly created domains. But these stable domain walls create problems in standard cosmology because if the energy density of the stable domain walls starts to dominate the energy density of the Universe, the Universe would have a rapid expansion (∝ t 2 ) which is incompatible with standard cosmology. Hence to avoid the domination of domain walls one can adopt a bias term [13] to the scalar potential to destabilise the domain walls and eventually collapse them. This term creates an energy difference between the degenerate vacua and this energy difference affects the domain walls with a volume pressure. When this volume pressure exceeds the surface tension of the wall, annihilation of the domain walls takes place. During the process of annihilation of domain walls, some parts of its energy are converted to GWs and contribute as a stochastic background of GW in the present Universe.

Domain Walls Formation
In this work we consider planar domain walls perpendicular to the z axis [90] to be formed around the boundaries of the domains after the spontaneous breaking of the discrete symmetries Z and Z . We define the scalars S 1 and S 2 as where we fix the radial part to their minima and introduce two phase factors φ 1 and φ 2 . The potential in Eq. 1 can be rewritten as The kinetic part of the Lagrangian can be expressed in terms of φ 1 and φ 2 as Using the potential in Eq. 31 the two field equations for φ 1 and φ 2 can be obtained as In order to obtain the domain wall solution we solve the above two equations (Eqs. 35 and 36) with the boundary conditions φ i → 2πn 2 at z → −∞ and φ i → 2π (n + 1) 2 at z → ∞ (with i = 1, 2 and n=0, 1, 2) [90]. By considering v 2 1 = v 2 2 (for analytical simplification) we get the following solution for n = 0, are equal. The width of the domain wall can be computed from the domain wall solution as The expression for the energy density ρ wall of the domain walls is given by Note that, we subtracted a constant term to satisfy the condition ρ wall → 0 for z → ±∞. The domain wall tension σ wall can be computed by integrating the energy density ρ wall over the z axis:

GWs Production from Domain Walls
As discussed earlier, the unstable domain walls would collapse and produce a significant amount of GWs. In this section, we briefly present the calculations to obtain the energy densities of GWs produced from the annihilation of domain walls [91,92].
It is discussed before that domain walls must collapse or annihilate before they can overclose the Universe. The annihilation occurs when the volume pressure force tends to overcome the tension force. From this condition one can calculate the annihilation time using the expression [13] t ann = 6.58 where A = 0.8 ± 0.1 is the area parameter [93] and C ann is a model depended parameter of order O(1). Here we take C ann = 5 for Z 2 symmetry [94]. In the above equation (Eq. 41), V bias denotes the bias term of the potential introduced to make the domain walls unstable. The temperature at the annihilation time t = t ann can be defined as [13] T ann = 3.41 where g * (T ann ) is the relativistic degrees of freedom for the energy density at temeperatute T ann . From Eq. 41, it is clear that t ann depends on the bias term. Hence to annihilate the domain walls before they can overclose the Universe (t ann < t dominate ), we require a lower bound on the magnitude of V bias [13] This condition also constrains the value of the annihilation temperature T ann [13] T ann > 1.62 If the domain walls decay into the SM particles, the decay products would affect the creation of light elements at the era of BBN, which is not favourable by the present observational results. Therefore it is required that the lifetime of domain walls should be shorter than t ann 0.1 sec [95,96]. This condition provides another bound on the magnitude of V bias and this bound is given by [13] The peak intensity of GWs (produced from the annihilations of domain walls) at the present time t 0 is given by [13] Ω dw GW h 2 (t 0 ) peak = 7.2 × 10 −18 where GW = 0.7 ± 0.4 [93] is the efficiency parameter and g * (T ann ) is the relativistic degrees of freedom at annihilation temperature for the entropy density. It can be mentioned here that in our analysis we assume that the GWs are produced from the annihilation of domain walls in the radiation dominated era. The peak frequency of the GWs for this scenario can be estimated from the following expression [13] f dw (t 0 ) peak = 1.1 × 10 −9 Hz g * (T ann ) 10 In Ref. [97], the authors mentioned that the domain walls cannot be formed if the bias term V bias is large enough. Hence an upper bound V bias << V should also be considered. We take V bias as the free parameter of the model. By varying V bias and by calculating model-dependent domain wall tension we compute the GW intensities. Needless to mention that when we choose V bias the constraints on V bias expressed in Eqs. 43, 45 are respected.
There are some recent studies on GW from domain walls based on numerical simulations that show Ω dw GW h 2 ∝ f 3 for f < f dw peak and Ω dw GW h 2 ∝ 1 f for f f dw peak [91]. In this work we use these relations and choose five sets of BPs (related to our model parameters) to calculate the GW intensity and peak frequency in the case of the present particle physics model.

Production from Strong First-Order Phase Transition
In this section, we pursue the electroweak phase transition and possible production mechanism of GWs from SFOPT in our proposed two complex scalars extended SM.

Finite Temperature Effective Potential
In order to explore the electroweak phase transition (EWPT) in the present model we include two potential terms with the tree-level potential V tree−level . Now the modified finite temperature effective potential can be expressed as [98] where V T =0 1−loop is the zero temperature Coleman-Weinberg one-loop effective potential and V T =0 1−loop is the finite temperature one-loop effective potential. The zero-temperature tree-level potential We obtain the above equation (Eq. 49) from Eq. 1 by replacing the scalar fields H, S 1 , S 2 with their VEVs v, v 1 , v 2 respectively. The zero-temperature Coleman-Weinberg one-loop effective potential can be written as [98] where the '+' and '−' symbols refer to the sign of the bosons and fermions. The summation i is over all the particles associated in the model and i ≡ (h 1 , h 2 , h 3 , G ± , G 0 , χ 1 , χ 2 , W, Z, t). In Eq. 50 n i , m i and C i denote the number of degrees of freedom, the field-dependent masses and renormalisation-scheme-dependent numerical constant of the particle i respectively. The degrees of freedom of the above mentioned particle species are (n W ± ) L = 4, (n W ± ) T = 2, (n Z ) L = 2, (n Z ) T = 1, n t = 12, n G ± = 2 and n h,h 2 ,h 3 ,G 0 ,χ 1 ,χ 2 = 1.
where g, g and y t refer to the SU (2) L gauge coupling, U (1) Y gauge coupling and top Yukawa coupling of the SM respectively. We perform our analysis by considering the Landau gauge where the Goldstone bosons are massless at zero temperature theory and also consider that ghost contributions are not appearing here [99]. In Eq. 50, the quantity Q is the renormalisation scale which we fix as Q = v = 246.22 GeV in our calculations. The finite temperature one-loop effective potential is given by [98] with Here we would like to mention that the values of the classical VEVs i.e (v, v 1 , v 2 ) change with temperature and at T = 0 it tends to the classical fixed values. We apply daisy resummation method [100] in our work for executing the thermal correction to the boson masses as µ 2 H (T ) = µ 2 H + c 1 T 2 , µ 2 S 1 (T ) = µ 2 S 1 + c 2 T 2 and µ 2 S 2 (T ) = µ 2 S 2 + c 3 T 2 , where and For the gauge bosons there are only finite temperature corrections to the longitudinal components and the corrected thermal masses of W and Z boson are given by [99] and We use the publicly available package CosmoTransitions [98] to include the finite temperature corrections with the tree-level potential.

Gravitational Waves Production from Strong First-Order Phase Transition
In this section, we discuss how the gravitational waves are produced from the SFOPT principle. Initially, the Universe is in a global minimum and with the decrease in temperature as the Universe evolves a minimum develops (metastable state). When the newly formed minima as a function of temperature tends to the true minimum then the process of tunnelling from false (global) to true minimum takes place through the nucleation of bubbles. During the process, a situation comes when the two minima become degenerate and that temperature indicates the critical temperature T c . At the nucleation temperature (less than T c ) the phase transition completes where at least one bubble is formed per unit time per Hubble volume. The latent heat energy liberated during the phase of first-order phase transition process contributes to the energy density of the stochastic GW background. The production of GWs is principally caused by three mechanisms namely (i) bubble collisions, (ii) sound waves induced by the bubbles running through the cosmic plasma and (iii) turbulence induced by the bubble expansions in the cosmic plasma.
The probability of formation of a bubble per unit time per unit volume at a temperature T can be written as [101] Γ = Γ 0 (T ) e −S 3 (T )/T , where Γ 0 (T ) scales as Γ 0 (T ) ∝ T 4 and the Euclidean action of the critical bubble S 3 (T ) is given by [101] S 3 = 4π dr r 2 1 2 where φ = (h, s 1 , s 2 ) represents a vector of the scalar fields in the potential V (Eq. 1) and V ef f is the finite temperature effective potential expressed in Eq. 48. Nucleation of bubbles occurs at a temperature T n where T n is the temperature where the conditions Γ ∼ 1 and S 3 (T n ) /T n ≈ 140 [98] are obeyed. As mentioned, the production mechanism of GWs from the strong first-order electroweak phase transition are driven by three processes namely bubble collisions [18]- [23], sound waves induced by the bubbles running through the cosmic plasma [24]- [27] and turbulence induced by the bubble expansions in the cosmic plasma [28]- [32]. The total GW intensity Ω GW h 2 from the SFOPT for a particular frequency f can be obtained by adding the contributions of the above three mechanisms and we have [18]- [33] The contribution of GW intensity from the bubbles collision Ω col h 2 is given by where β is the inverse time scale of the phase transition parameter and it has the following form where H denotes the Hubble parameter at the nucleation temperature T n . We estimate the bubble wall velocity v w using the following expression as [102,22,48,49] In Eq. 64 the quantity κ is the fraction of latent heat energy deposited in a thin shell which takes the form with [52,103] α ∞ = 30 where v n is the VEV of Higgs at the nucleation temperature T n , m W , m Z and m t denotes the masses of the gauge bosons W , Z and the top quark t respectively. The parameter α can be defined as the ratio of energy density difference between false and true vacuum released during the electroweak phase transition ρ vac to the background energy density of the plasma ρ rad * at T n . The expression of α is given by with and The expression of peak frequency f col (Eq. 64) produced by the bubble collisions is the suppression factor turns out to be less than 1, then we need to include the factor to the sound wave component of the GW intensity. The contribution Ω turb h 2 from the turbulence in the plasma is given by where = 0.1 and the peak frequency f turb for the turbulence mechanism reads: where the parameter h * = 16.5 × 10 −6 Hz T n 100 GeV We use Eqs. 63 -78 for the estimation of the gravitational wave intensity in this work.

Calculations and Results
In this section, we estimate the intensities of the gravitational wave produced via domain wall annihilation as well as from strong first-order phase transition for the present extended SM. To check its detectability, our calculated model-dependent GW intensities are compared with the future space-based and ground-based interferometers such as ALIA, BBO, DECIGO, TianQin, Taiji, aLIGO, aLIGO+ and pulsar timing arrays such as SKA, IPTA, EPTA, PPTA, NANOGrav11 and NANOGrav12.5. For the calculation of GW intensities we choose five benchmark points (BPs) such that the numerical values of the model parameters in each of the chosen BPs satisfy the theoretical constraints (such as vacuum stability, perturbativity) as well as the experimental constraints (such as LHC). In case of GW production from domain walls, GW intensity depends mainly on factors such as domain wall tension σ wall , annihilation time t ann , annihilation temperature T ann , degrees of freedom at annihilation time g * s (T ann ) and bias term V bias . On the other hand for the GW production via SFOPT, the strength of the first-order phase transition parameter α, inverse of the time-scale of the phase transition parameter β, bubble wall velocity v w , nucleation temperature T n and value of Higgs VEV v n at T n play significant roles. Using the equations given in section 4, we compute and explore the formation of domain walls and production of GW from the annihilation of the domain walls for the present model considered in this work. In order to calculate the domain wall tension we first solve the Eqs. 35 and 36 and get the domain wall solution. Then we estimate the domain wall tension using Eq. 39 and 40. The domain walls are then made unstable by introducing a bias term V bias . Here we take V bias as a free parameter. We also use two bounds on V bias to make the domain walls annihilate before it could overclose the Universe and before BBN era. Considering these bounds we compute the annihilation time of the domain walls and annihilation temperature at that time by using Eq. 41 and Eq. 42 respectively. Finally, we calculate the GW intensity and the peak frequency at the present time from the annihilation of domain walls using Eqs. 46 and 47 respectively. Note that, our chosen five BPs for the calculation of GW intensity satisfy the bounds on V bias and T ann .  Table 1: The chosen five benchmarks points (BPs, BP 1-5) to explore the GW production from both the domain walls and strong first-order phase transition in two complex scalars extended SM.
On the other hand, for the calculation of GW intensity from the first-order phase transition we consider a finite temperature effective potential (discussed in section 5). The finite temperature effective potential is obtained by adding two terms with the tree-level potential -one is the oneloop CW potential at zero temperature V T =0 1−loop and the other is the finite temperature potential V T =0 1−loop . We consider the three possible processes such as bubble collisions, sound waves and turbulence in the plasma for the production of GWs from phase transition. Then we compute the total GW intensity using Eqs. 63 -78 from the SFOPT. We also mention that in this work we choose five BPs to demonstrate that our model can induce strong first-order phase transition for each of the five choices which is a necessary condition for the production of GWs [99,106]. The choice of the five BPs with the values of the model parameters are furnished in Table 1.
Here we use a package namely Cosmotransition for calculating the thermal parameters related to the phase transition properties for all the chosen BPs.
We tabulate our calculated results in Tables 2-5. While the calculations related to the GWs from domain walls are furnished in Table 2 and Table 3, in Table 4 and Table 5 we show the computed results when GWs are formed via SFOPT. In Table 2 we show our results for the quantities g * (T ann ), g * s (T ann ), t ann , T ann , V bias , σ wall (related to the estimation of GW intensity from domain walls) for the five chosen BPs. We present our calculated GW intensity at the present epoch and the corresponding peak frequency for all the BPs (BP 1-5) in Table 3.
The values of the thermal parameters (v c , T c , v n , T n , α, β/H) which are the useful parameters for the calculation of GW intensity from SFOPT are tabulated in Table 4 for all the five BPs. In Table 4 we also indicate the critical temperature T c (the temperature at which two phases are degenerate) and the numerical value of the v c , the VEV at that temperature T c . Note that all the BPs satisfy the condition v c /T c 1 except the ones in BP 4. For BP 4 v c /T c < 1, so in that case the SFOPT condition is not satisfied and hence it is unable to produce GWs from SFOPT. So in further calculations of GW intensity we do not calculate the intensity for BP 4. Table  4 shows that the nucleation temperature T n is less than T c . Therefore with further decrease in temperature from T c , nucleation of the bubble occurs and phase transition is completed. In Table   4, we also present the value of HR * U f , the factor which indicates whether or not the contribution of sound wave is survives more than Hubble time. In case HR * U f < 1, the SW contribution decays within the Hubble time. Following Refs. [103,104,105] we compute this factor for all BPs and obtain HR * U f < 1 for all the BPs (BP 1, BP 2, BP 3, BP 5) except BP 4. Therefore as mentioned in Section 5 we multiply the factor with the sound wave component of GW intensity (Eq. 73). In Table 5 we show the results for GW intensity and corresponding peak frequency from SFOPT. We obtain a single peak in GW intensity for the case of each BP but when BP 3 is adopted, two peaks appears instead for GW intensity when varied with GW frequency. The reason behind the two peaks being two GW components, bubble collisions and SW, play dominant roles. While in the other cases only bubble collision plays a major role and only bubble collision dominantly influences the peak of the GW spectrum. In case of BP 1, 2 and 5 the peak frequencies of the total GW spectrum from SFOPT are equal to the peak frequencies of the component of bubble collision f col = f P T peak . However in the case of BP 3 for which two peaks are obtained, one peak corresponds to f col and while the other is at f SW . Note that, since BP 4 does not satisfy the SFOPT criteria, it can not produce GWs from strong first-order phase transition.
C ann =5, A=0.8, GW =0.7. BP g * (T ann ) g * s (T ann ) t ann T ann V bias σ wall in sec in GeV in MeV 4 in TeV 3     The calculated values of peak frequency and its corresponding GW intensity from the contribution of strong first-order phase transition for the chosen BPs (BP 1, BP 2, BP 3 and BP 5). We obtain two peak frequencies for BP 3 and its corresponding GW intensities are also mentioned in Table 5. Note that we do not show the results for BP 4 since in the case of BP 4 the strong first-order phase transition condition is not satisfied (v c /T c < 1, Table 4), as a result it is unable to produce GWs.
In Figure 1, we show the variation of GW intensity from domain walls as a function of V bias with a fixed value for domain wall tension σ wall . In Figure 1 we adopt the benchmark values in BP 3 to demonstrate the variation of GW intensity with bias potential V bias . From Figure 1 one observes that smaller the value of V bias , more is the strength of GW intensity. That smaller values of V bias lead to larger intensity of GW can also be concluded from Eqs. 42 and 46. We repeat the computation of the variations of the GW intensity with V bias for all the other BPs and obtain similar results. Therefore we choose small values of V bias which satisfy the bounds mentioned in section 4.2, to calculate the GW intensity in case of GW production from domain walls. The selected values of V bias for all the BPs are mentioned in Table 2. In Figure 2 we show the variation of the parameter S 3 /T with temperature T for BP 3. We obtained similar variations for other BPs but for the purpose of demonstration we only present here the results for BP 3. By calculating the slope of the (S 3 /T vs T ) graph we compute the value of the parameter β using Eq. 65. In Figure 2 we also draw a horizontal line (red color) at S 3 /T = 140 to indicate the onset of bubble nucleation.
Finally in Figure 3 we plot the variation of GW intensity as a function of frequency. In the left panel of Figure 3 we show the GW intensity from domain walls for our selected five BPs whereas in the right panel of Figure 3 the GW intensity from phase transition for the same BPs are plotted. We perform a direct comparison between our calculated model-dependent GW intensity with the sensitivity plots [107] of PTAs such as SKA, IPTA, EPTA, PPTA, NANOGrav11 and NANOGrav12.5 for the domain wall case and found that the calculations using BP 4 lie within the sensitivity curves of SKA and IPTA. We note that calculations using BP 4 yield higher GW intensity compared to the results when other BPs are used when GWs from domain walls are considered but the GWs from first-order phase transition cannot be obtained if BP 4 is used for the calculation. We also find that the VEVs v 1 , v 2 play important roles for GW intensities from domain walls. From Tables 1, 3 and Figure 3 we may conclude that higher the values of VEVs, higher are the GW intensities for domain wall scenario.
We also compare our calculated GW spectrum from SFOPT (for same BPs) with the future space-based and ground-based gravitational wave detectors such as ALIA, BBO, DECIGO, LISA, TianQin, Taiji, aLIGO and aLIGO+ and find that GW intensities obtained using BP 1 lies within the sensitivity curves 1 of ALIA, BBO, DECIGO, Taiji and marginally lies in TianQin while for the case of BP 2 the GW intensities lie within the sensitivity curves of ALIA, BBO, DECIGO, Taiji. The GW intensities marginally lie within the sensitivity curves of ALIA, BBO, DECIGO when BP 3 is used and for BP 5 they lie within the sensitivity curves of ALIA, BBO and DECIGO.
From the results mentioned in Table 2 -5, we conclude that the dominant parameter influencing GW in case of GWs from domain wall is σ wall whereas in the case of SFOPT the dominant parameter is β.

Summary and Conclusions
In this work, we have explored two possible production mechanisms of GWs, one is from the annihilation of domain walls and the other is from the strong first-order phase transition in the early Universe, within the framework of a particle physics model where SM is extended by two complex scalars. We investigate the detection possibilities of these GWs at the pulsar timing arrays such as SKA, IPTA, EPTA, PPTA, NANOGrav11 and NANOGrav12.5 as well as at future space-based and ground-based gravitational wave detectors such as ALIA, BBO, DECIGO, LISA, TianQin, Taiji, aLIGO and aLIGO+. A discrete Z 2 × Z 2 symmetry is imposed in the model and the two added complex scalar fields S 1 and S 2 acquire a non zero vacuum expectation value when this imposed symmetry is spontaneously broken. As the discrete symmetry is broken spontaneously, domain walls are formed which is made unstable by considering this symmetry to be approximate which is explicitly broken by introducing a bias term in the theory. After the formation of these unstable domain walls, they annihilate and produce a significant amount of GWs. We also discuss the production procedure of GWs from a strong first-order phase transition. For that, we consider a finite temperature effective potential and show that the potential induces a first-order phase transition. We choose five BPs to calculate the GW intensity and frequency from both the production mechanisms (annihilation of domain walls and strong first-order phase transition). The BPs are chosen in such a way that they satisfy both the theoretical constraints such as vacuum stability, perturbativity and the experimental constraints such as collider bounds. We found that the VEVs v 1 , v 2 (VEVs acquired by S 1 and S 2 after spontaneous breaking of Z and Z respectively) play a very important role to calculate the GW intensity from the collapse of domain walls and higher values of VEVs give higher GW intensities. The peaks for the GW intensities for their production from two mechanisms, namely domain wall annihilation and strong first-order phase transition, appear in two different frequency regions. For the case of GWs from domain wall annihilation, the intensity peaks in a lower frequency regime around ∼ 10 −9 Hz whereas for GWs from strong first-order phase transition, these peaks are obtained at comparatively higher frequency regime of ∼ 10 −3 −10 −2 Hz. We also demonstrate that our calculated GWs from the annihilations of domain walls can be probed by low-frequency PTAs such as SKA and IPTA whereas the GWs which are produced from strong first-order phase transition can be accessed by comparably higher frequency detectors such as ALIA, BBO, DECIGO, TianQin and Taiji. In this work, within the framework of our proposed model we show that the simple extension of SM can explain the formation of unstable domain walls as well as the strong first-order electroweak phase transition. One can also extend this model for a particle dark matter theory by considering one of the components of the complex scalars to be a dark matter candidate and work out its phenomenology as well as the viability of such a dark matter candidate. This is for posterity.