# Left-Right Symmetry Breaking and Gravitational Waves : A Tale of Two Phase
Transitions

Z. A. Borboruah  U. A. Yajnik

###### Abstract

We study possible ways gravitational waves (GW) are sourced in a theory with
minimal left-right symmetry breaking. Generically first order phase
transitions (FOPT) lead to gravitational waves sourced by bubble collisions,
while second order phase transitions (SOPT) do not. Interesting variants on
the standard classification of phase transitions occur due to the breaking of
discrete parity combined with the limitation of light cone in the early
Universe. If local effective potential signals SOPT or a cross over, breaking
of discrete parity in conjunction with the finiteness of the causal horizon
leads to a causal horizon limited second order phase transition, which results
in domain walls (DWs) separating left-like and right-like domains. Thus the
putative SOPT also gives rise to GW via decaying domain walls, testable at
experiments such as IPTA, DECIGO, and LISA. On the other hand, for the case of
FOPT, we get the usual signal from spontaneously created bubbles, but also a
late forming domain wall structure separating the two types of vacua. Thus the
traditional FOPT case gives rise to two distinct peaks in the spectrum of GW,
verifiable for the low symmetry breaking scales $10^{4}-10^{6}$ GeV, but
beyond the reach of currently planned experiments for a high scale $\sim
10^{10}$GeV. Finally, we point out that a version of the left-right symmetric
model which separates parity breaking from gauge symmetry breaking is also
subject to domain wall formation and amenable to GW observations.

##  1 Introduction

The detection of gravitational waves (GWs) [1] has opened a new window into
the unknown of the Universe. Originating from the linearized theory of general
relativity [2], GWs are a new independent source of information about cosmic
events. Also, since gravitation is the weakest force in nature, gravitons are
believed to decouple from the rest of the forces and matter as early as the
Planck scale itself, meaning they can convey information from the very dawn of
the Universe.

The early Universe GW sources include not only astronomical objects, but
cosmological events like inflation, phase transitions, etc. Throughout time,
GWs from all such events superimpose to create what’s called a _stochastic
background of GWs_ , which is believed to be homogeneous, isotropic,
unpolarized with an amplitude following a Gaussian distribution [3]. Several
current and future experiments are aimed at detecting this background,
directly via interferometers [4, 5] or through indirect methods like Pulsar
Timing Arrays (PTAs) [6]. This GW background can be a good alternative probe
for high and ultra high energy phenomena including those arising from grand
unified theories (GUTs) and other particle physics theories beyond the
Standard Model (SM).

One such high energy theory is the Left-Right Symmetric Model (LRSM),
especially since it accommodates massive neutrinos with well defined gauge
charges for the new right handed neutrino states. Based on the gauge group
$SU(2)_{L}\otimes SU(2)_{R}\otimes U(1)_{B-L}$, this model is fully symmetric
between right and left chiral states at high temperatures in the early
Universe when the thermal ground state does not break symmetry spontaneously.
As such why low energy physics is governed by $SU(2)_{L}\times U(1)_{Y}$ is
unknown. An early proposal to rectify this situation was made in [7, 8] by
dissociating spontaneous breaking of the so called _D-parity_ from the
breaking of right-handed gauge symmetry, by invoking an extra scalar beyond
the basic LRSM particle content.

But as we will discuss later in this paper, this scenario also cannot avoid
domain walls and will also contribute to the GW background. Conventionally
LRSM is spontaneously broken down to SM through additional Higgs particles,
for example using scalar $SU(2)$ doublets or triplets in the supersymmetric
(SUSY) [9, 10] and Non-SUSY [11, 12, 13, 14, 15, 16, 17, 18, 19] contexts. A
minimal left-right symmetric model (MLRSM) with scalar triplets has garnered
much attention because it can implement the seesaw mechanism and can be
embedded in a $SO(10)$ grand unified theory (GUT) nicely [20, 21, 22, 23, 24].
LRSM has been studied extensively in view of the present and future colliders
[25, 26, 27, 28, 29, 30, 31, 32] but LHC has detected no evidence of it up to
$\mathcal{O}(10\text{ TeV})$. Therefore the GW background presents the
exciting possibility of providing the first hints of this theory, especially
if originating at significantly higher energies, thus complementing the role
of future colliders.

The breaking of discrete symmetry $L\leftrightarrow R$ which is generically
exact in left-right symmetric models leads to the formation of domain walls
(DW) in the early Universe. Domain walls are massive topological defects
created when a discrete symmetry is broken [33, 34]. They separate two regions
of distinct degenerate vacua of the theory. This is a problem for standard
cosmology, as stable domain walls come to dominate the energy density of the
Universe at later times [35] and conflict with structure formation and Big
Bang Nucleosynthesis (BBN). Therefore, they must be transient [36, 37, 38, 39,
40, 41]. This can be achieved by making one of the vacua slightly lower in
energy so that the domain walls move under the vacuum pressure. Eventually,
violent collisions between DW cause their annihilation returning the Universe
to a homogeneous state after entropic equilibrium is reached.

While it is usual to identify phase transitions in gauge theories as First
Order (FOPT), Second order (SOPT), or crossover, the left-right symmetric
models present a peculiar situation, due to the presence of causal horizons in
the Big Bang Universe. While the effective potential may signal a second-order
phase transition, the correlation lengths remain limited by the causal
horizon, thus evading the divergence of length scales characteristic of the
SOPT. The situation is further complicated by the fact that there are two
distinct fields whose vacuum expectation values (VEVs) signal degenerate
vacua. In causally disconnected regions distinct fields acquire VEV on either
side of the wall, eg., $SU(2)_{L}$ triplet on one side of the wall (right-like
vacuum), and $SU(2)_{R}$ triplet (left-like vacuum) on the other. This is also
the case when the local effective potential signals a cross over phase
transition. This situation is best characterized as “causal horizon limited
SOPT”, to be abbreviated as CHL-SOPT. The eventual state resulting from such a
phase transition will be a network of domain walls, with the domains
alternately enclosing left-like and right-like vacua.

Likewise, even the putative FOPT is fundamentally modified in that the
spontaneous formation of bubbles of true vacuum could either lead to bubbles
of left-like vacuum or right-like vacuum. When these bubbles grow and
encounter each other, the collision can result in a larger bubble, provided
the interiors contain the same type of vacuum. On the other hand, if one
bubble contains left-like while the other contains right-like vacuum then the
result is a stable left-right domain wall. The wall configuration is made up
of non-trivial condensates of two or more fields [42, 43]. The situation is
depicted in Fig. 1. The eventual result of such a phase transition will also
be a network of domain walls enclosing distinct types of vacua, making it look
similar to SOPT of a different scale. We refer to this as the Degenerate Field
FOPT.

Gravitational waves from first-order phase transition [44, 45, 46, 47, 48, 49,
50] and more specifically from MLRSM [51, 52] have been studied extensively
which we will discuss in section 5. The several contributions to the GW
resulting from SOPT have been identified to be [53, 54, 55] from 1) evolution
of the initial quantum fluctuations of the fields and 2) late time evolution
of massive Domain Walls. At late times, the contribution from the former
becomes insignificant compared to the contribution from the latter. This paper
will consider the CHL-SOPT case, and calculate the energy spectrum of the
gravitational waves produced. Our approach parallels the methods used in [53,
54, 55] but with the difference that we consider the late time evolution of
domain walls after the Ginzburg temperature has been reached in the context of
LRSM. For the Degenerate Field FOPT case, we determine the double peak
spectrum for models with symmetry breaking scales $10^{4}-10^{6}$ GeV and
$\sim 10^{10}$GeV and discuss the possibility of their verifiability.

The arrangement of this paper is as follows: in the next section, Sec. 2 we
describe the MLRSM with triplets and the loop corrected temperature dependent
effective potential. In a subsection we emphasize the cosmology of the two
kinds of phase transitions, identifying CHL-SOPT and Degenerate Field FOPT. In
Sec. 3, we discuss GWs from CHL-SOPT followed by Sec. 4 describing the
numerical methods used to simulate the decaying domain walls. In Sec. 5, we
discuss GWs for the Degenerate Field FOPT case, identifying the two peak
feature. In sec. 6 we show that the spontaneous parity breaking scenario of
[8, 7] also results in characteristic GW background. Sec. 7 contains a summary
and conclusion. Appendix A contains some of the additional details of the
numerical procedure used.

##  2 Minimal Left-Right Symmetric Model

The MLRSM [56, 57, 58, 59, 60, 61, 62, 63] is an appealing extension of the SM
which naturally accords gauge charges to the right-handed neutrino states
which are generically needed to provide masses to neutrinos. It is based on
the gauge group,

$\displaystyle\mathcal{G}_{LR}\equiv SU(2)_{L}\times SU(2)_{R}\times
U(1)_{B-L}\times SU(3)_{C}\,$ (2.1)

where the electric charge is defined as

$\displaystyle Q=T_{3L}+T_{3R}+\frac{B-L}{2}\,.$ (2.2)

Under this left-right symmetric gauge group, the usual quarks and leptons
transform as

$\displaystyle q_{L}=\begin{pmatrix}u_{L}\\\
d_{L}\end{pmatrix}\equiv[2,1,1/3,3]\,,~{}q_{R}=\begin{pmatrix}u_{R}\\\
d_{R}\end{pmatrix}\equiv[1,2,1/3,3]\,,$
$\displaystyle\ell_{L}=\begin{pmatrix}\nu_{L}\\\
e_{L}\end{pmatrix}\equiv[2,1,-1,1]\,,~{}\ell_{R}=\begin{pmatrix}\nu_{R}\\\
e_{R}\end{pmatrix}\equiv[1,2,-1,1]\,.$

The Higgs sector may vary in LRSM. We will focus on the minimal Higgs sector
containing the triplets and the bidoublet [64], $\Delta_{L}(1,0,2,0)$,
$\Delta_{R}(0,1,2,0)$, $\phi(1/2,1/2^{*},0,0)$. They can be represented in the
matrix form,

$\phi=\left(\begin{array}[]{cc}\phi_{1}^{0}&\phi_{1}^{+}\\\
\phi_{2}^{-}&\phi_{2}^{0}\end{array}\right),\quad\Delta_{L}=\left(\begin{array}[]{cc}\frac{1}{\sqrt{2}}\delta_{L}^{+}&\delta_{L}^{++}\\\
\delta_{L}^{0}&-\frac{1}{\sqrt{2}}\delta_{L}^{+}\end{array}\right),\quad\Delta_{R}=\left(\begin{array}[]{cc}\frac{1}{\sqrt{2}}\delta_{R}^{+}&\delta_{R}^{++}\\\
\delta_{R}^{0}&-\frac{1}{\sqrt{2}}\delta_{R}^{+}\end{array}\right)$ (2.3)

so that they transform as

$\displaystyle SU(2)_{L}\times SU(2)_{R}$ $\displaystyle:\phi\rightarrow
U_{L}\phi U_{R}^{\dagger},$ $\displaystyle\Delta_{L}\rightarrow
U_{L}\Delta_{L}U_{L}^{\dagger},\quad\Delta_{R}\rightarrow
U_{R}\Delta_{R}U_{R}^{\dagger}$ (2.4) $\displaystyle U(1)_{B-L}$
$\displaystyle:\phi\rightarrow\phi,$ $\displaystyle\Delta_{L}\rightarrow
e^{i2\theta}\Delta_{L},\quad\Delta_{R}\rightarrow e^{i2\theta}\Delta_{R}$
(2.5)

for general symmetry transformations: $U_{L}\in SU(2)_{L},U_{R}\in SU(2)_{R}$,
and $e^{i\theta}\in U(1)_{B-L}$.

###  2.1 Potential

The tree-level standard scalar potential of the LRSM with triplets is given
by,

$V_{0}=V_{\phi}+V_{\Delta}+V_{\phi\Delta}$ (2.6)

where

$\displaystyle V_{\phi}=$
$\displaystyle-\mu_{1}^{2}\operatorname{Tr}\left[\phi^{\dagger}\phi\right]-\mu_{2}^{2}\left(\operatorname{Tr}\left[\tilde{\phi}\phi^{\dagger}\right]+\operatorname{Tr}\left[\tilde{\phi}^{\dagger}\phi\right]\right)+\lambda_{1}\operatorname{Tr}\left[\phi^{\dagger}\phi\right]^{2}+\lambda_{2}\left(\operatorname{Tr}\left[\tilde{\phi}\phi^{\dagger}\right]^{2}+\operatorname{Tr}\left[\tilde{\phi}^{\dagger}\phi\right]^{2}\right)$
$\displaystyle+\lambda_{3}\operatorname{Tr}\left[\tilde{\phi}\phi^{\dagger}\right]\operatorname{Tr}\left[\tilde{\phi}^{\dagger}\phi\right]+\lambda_{4}\operatorname{Tr}\left[\phi^{\dagger}\phi\right]\left(\operatorname{Tr}\left[\tilde{\phi}\phi^{\dagger}\right]+\operatorname{Tr}\left[\tilde{\phi}^{\dagger}\phi\right]\right),$
$\displaystyle V_{\Delta}=$
$\displaystyle-\mu_{3}^{2}\left(\operatorname{Tr}\left[\Delta_{L}\Delta_{L}^{\dagger}\right]+\operatorname{Tr}\left[\Delta_{R}\Delta_{R}^{\dagger}\right]\right)+\rho_{1}\left(\operatorname{Tr}\left[\Delta_{L}\Delta_{L}^{\dagger}\right]^{2}+\operatorname{Tr}\left[\Delta_{R}\Delta_{R}^{\dagger}\right]^{2}\right)$
$\displaystyle+\rho_{2}\left(\operatorname{Tr}\left[\Delta_{L}\Delta_{L}\right]\operatorname{Tr}\left[\Delta_{L}^{\dagger}\Delta_{L}^{\dagger}\right]+\operatorname{Tr}\left[\Delta_{R}\Delta_{R}\right]\operatorname{Tr}\left[\Delta_{R}^{\dagger}\Delta_{R}^{\dagger}\right]\right)+\rho_{3}\operatorname{Tr}\left[\Delta_{L}\Delta_{L}^{\dagger}\right]\operatorname{Tr}\left[\Delta_{R}\Delta_{R}^{\dagger}\right]$
$\displaystyle+\rho_{4}\left(\operatorname{Tr}\left[\Delta_{L}\Delta_{L}\right]\operatorname{Tr}\left[\Delta_{R}^{\dagger}\Delta_{R}^{\dagger}\right]+\operatorname{Tr}\left[\Delta_{L}^{\dagger}\Delta_{L}^{\dagger}\right]\operatorname{Tr}\left[\Delta_{R}\Delta_{R}\right]\right),$
$\displaystyle V_{\phi\Delta}=$
$\displaystyle\alpha_{1}\operatorname{Tr}\left[\phi^{\dagger}\phi\right]\left(\operatorname{Tr}\left[\Delta_{L}\Delta_{L}^{\dagger}\right]+\operatorname{Tr}\left[\Delta_{R}\Delta_{R}^{\dagger}\right]\right)+\alpha_{3}\left(\operatorname{Tr}\left[\phi\phi^{\dagger}\Delta_{L}\Delta_{L}^{\dagger}\right]+\operatorname{Tr}\left[\phi^{\dagger}\phi\Delta_{R}\Delta_{R}^{\dagger}\right]\right)$
$\displaystyle+\alpha_{2}\left(\operatorname{Tr}\left[\Delta_{L}\Delta_{L}^{\dagger}\right]\operatorname{Tr}\left[\tilde{\phi}\phi^{\dagger}\right]+\operatorname{Tr}\left[\Delta_{R}\Delta_{R}^{\dagger}\right]\operatorname{Tr}\left[\tilde{\phi}^{\dagger}\phi\right]+\text{
h.c. }\right)$
$\displaystyle+\beta_{1}\left(\operatorname{Tr}\left[\phi\Delta_{R}\phi^{\dagger}\Delta_{L}^{\dagger}\right]+\operatorname{Tr}\left[\phi^{\dagger}\Delta_{L}\phi\Delta_{R}^{\dagger}\right]\right)+\beta_{2}\left(\operatorname{Tr}\left[\tilde{\phi}\Delta_{R}\phi^{\dagger}\Delta_{L}^{\dagger}\right]+\operatorname{Tr}\left[\tilde{\phi}^{\dagger}\Delta_{L}\phi\Delta_{R}^{\dagger}\right]\right)$
$\displaystyle+\beta_{3}\left(\operatorname{Tr}\left[\phi\Delta_{R}\tilde{\phi}^{\dagger}\Delta_{L}^{\dagger}\right]+\operatorname{Tr}\left[\phi^{\dagger}\Delta_{L}\tilde{\phi}\Delta_{R}^{\dagger}\right]\right)$
(2.7)

where $\tilde{\phi}=\sigma_{2}\phi^{*}\sigma_{2}$. All the couplings are
assumed to be real for simplicity. The vacuum expectation values (VEVs) of the
Higgs fields that give the correct breakdown of LR symmetry are as follows,

$\langle\phi\rangle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}\kappa_{1}&0\\\
0&\kappa_{2}e^{i\theta_{2}}\end{array}\right),\quad\left\langle\Delta_{L}\right\rangle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}0&0\\\
v_{L}e^{i\theta_{L}}&0\end{array}\right),\quad\left\langle\Delta_{R}\right\rangle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}0&0\\\
v_{R}&0\end{array}\right)$ (2.8)

The bidoublet VEVs provide the electroweak scale,

$\sqrt{\kappa_{1}^{2}+\kappa_{2}^{2}}=v=246\text{ GeV}$ (2.9)

and the parameter $\beta=\arctan{\kappa_{2}/\kappa_{1}}$ captures the ratio of
the VEV’s. The minimization conditions on the potential, with respect to the
four VEVs lead to a relationship between $v_{L}$ and $v_{R}$ through a see-
saw-like formula,

$\beta_{1}\kappa_{1}\kappa_{2}\cos\left(\theta_{2}-\theta_{L}\right)+\beta_{2}\kappa_{1}^{2}\cos\theta_{L}+\beta_{3}\kappa_{2}^{2}\cos\left(2\theta_{2}-\theta_{L}\right)=\left(2\rho_{1}-\rho_{3}\right)v_{L}v_{R}$
(2.10)

It is customary to set $\beta_{1,2,3}=0$ so that the strong condition of this
see-saw can be circumvented, $v_{L}$ can be chosen $0$ while $v_{R}$ can be
assigned phenomenologically desirable large value. Additionally $\theta_{2}$
and $\theta_{L}$ may be set to $0$ for simplicity.

We shall use the VEVs from SM and expected phenomenology as the inputs in the
minimization conditions, in which case the parameters $\mu^{2}_{1,2,3}$ get
determined as follows.

$\begin{array}[]{l}\mu_{1}^{2}=\lambda_{1}\left(\kappa_{1}^{2}+\kappa_{2}^{2}\right)+2\kappa_{1}\kappa_{2}\lambda_{4}+\frac{1}{2}v_{R}^{2}\alpha_{1}-\frac{\alpha_{3}}{2}\frac{v_{R}^{2}\kappa_{2}^{2}}{\kappa_{1}^{2}-\kappa_{2}^{2}}\\\
\mu_{2}^{2}=\left(2\lambda_{2}+\lambda_{3}\right)\kappa_{1}\kappa_{2}+\frac{\lambda_{4}}{2}\left(\kappa_{1}^{2}+\kappa_{2}^{2}\right)+\frac{\alpha_{2}}{2}v_{R}^{2}+\frac{\alpha_{3}}{4}\frac{v_{R}^{2}\kappa_{1}\kappa_{2}}{\kappa_{1}^{2}-\kappa_{2}^{2}}\\\
\mu_{3}^{2}=\rho_{1}v_{R}^{2}+\frac{1}{2}\alpha_{1}\left(\kappa_{1}^{2}+\kappa_{2}^{2}\right)+2\alpha_{2}\kappa_{1}\kappa_{2}+\frac{1}{2}\alpha_{3}\kappa_{2}^{2}\end{array}$
(2.11)

An important point that has not received due attention in the literature is
that the see-saw-like relation equally well implies a possible large value for
$v_{L}$ and $v_{R}\approx 0$ due to the left-right symmetry of the model.
While low energy phenomenology has preferred the $v_{L}\approx 0$ solution,
this will not persist at high temperatures in the early Universe. Indeed we
expect gauge symmetry restoration at high temperature [65, 66, 67], a scale
$T_{c}$ above which the thermal ground state has no non-zero VEVs. When the
Universe cools it must choose between the left-like vacuum with $v_{L}\neq 0$
or the right-like vacuum with $v_{R}\neq 0$. As already introduced in Sec. 1,
this leads to intricate possibilities for the nature of the phase transition,
as we further discuss in Sec. 2.3

###  2.2 Effective Potential

From low energy phenomenology we are interested in scenarios where the
$SU(2)_{R}$ breaking scale is much higher than electroweak symmetry breaking
$v_{R}\gg\kappa_{1},\kappa_{2}$. It turns out that it is sufficient to focus
on large $\delta_{R}^{0}$ values in field space, compared to the remaining
neutral fields $\phi_{1}^{0},\phi_{2}^{0}$ and $\delta_{L}^{0}$ [51]. However,
in the context of the early Universe, the alternative possibility of $v_{L}$
being large while ignoring $v_{R}$ also arises. Hence we retain the two
contributions in the tree level potential from $\delta_{R}^{0}$ and
$\delta_{L}^{0}$

$V_{0}^{total}(r,l)=V_{0}(r)+V_{0}(l)$ (2.12)

where $r=\delta_{R}^{0}/\sqrt{2}$ and $l=\delta_{L}^{0}/\sqrt{2}$, and $V_{0}$
are both identical functional forms of their arguments,

$V_{0}(r)\approx\frac{\rho_{1}}{4}\left(r^{2}-\eta^{2}\right)^{2}$ (2.13)

where $\eta=\sqrt{\mu_{3}^{2}/\rho_{1}}$. Note that it is sufficient to
consider only the first quadrant $r>0,l>0$ of the $r-l$ plane. Indeed due to
the see-saw type condition, where $r\rightarrow 0$, $l\rightarrow v_{L}\neq 0$
and vice versa. As such we expect domain walls, separating regions of left-
like vacuum $l\neq 0,r=0$ and $l=0,r\neq 0$ with the wall constituting a
transition region where one VEV goes to zero and the other turns on.

Let us focus on generalizing the functional form (2.13) of the potential to
the case of non-zero temperature. The temperature-dependent effective
potential is given by [51],

$V_{\mathrm{eff}}(r,T)=V_{0}(r)+V_{\mathrm{CW}}(r)+V_{\mathrm{FT}}(r,T)+V_{\mathrm{D}}(r,T)$
(2.14)

where

$\displaystyle V_{\mathrm{CW}}(r)=\frac{1}{64\pi^{2}}$
$\displaystyle{\left[\sum_{i}m_{i}^{4}(r)\left(\log\frac{m_{i}^{2}(r)}{\mu^{2}}-\frac{3}{2}\right)+6m_{W_{R}}^{4}(r)\left(\log\frac{m_{W_{R}}^{2}(r)}{\mu^{2}}-\frac{5}{6}\right)\right.}$
$\displaystyle\left.+3m_{Z_{R}}^{4}(r)\left(\log\frac{m_{Z_{R}}^{2}(r)}{\mu^{2}}-\frac{5}{6}\right)-6m_{\nu_{R}}^{4}(r)\left(\log\frac{m_{\nu_{R}}^{2}(r)}{\mu^{2}}-\frac{3}{2}\right)\right]$
(2.15)

is the Coleman-Weinberg correction to the tree potential and,

$V_{\mathrm{FT}}(r,T)=\frac{T^{4}}{2\pi^{2}}\left[\sum_{i}J_{\mathrm{B}}\left(\frac{m_{i}^{2}(r)}{T^{2}}\right)+6J_{\mathrm{B}}\left(\frac{m_{W_{R}}^{2}(r)}{T^{2}}\right)+3J_{\mathrm{B}}\left(\frac{m_{Z_{R}}^{2}(r)}{T^{2}}\right)-6J_{\mathrm{F}}\left(\frac{m_{\nu_{R}}^{2}(r)}{T^{2}}\right)\right]$
(2.16)

is the temperature correction from finite temperature field theory, where
$J_{B/F}$ are integrals for bosons/fermions given by,

$J_{\mathrm{B}/\mathrm{F}}\left(r^{2}\right)=\int_{0}^{\infty}\mathrm{d}xx^{2}\log\left(1\mp\mathrm{e}^{-\sqrt{x^{2}+r^{2}}}\right)$
(2.17)

The mass terms $m_{i},m_{W_{R}},m_{Z_{R}},m_{\nu_{R}}$ are the field ($r$)
dependent scalar, gauge and fermion masses [51]. The $\mu$ in (2.2) is the
renormalization scale which we can conveniently assume to be equal to $v_{R}$.
The last term includes the Daisy diagrams for higher loop corrections,

$V_{\mathrm{D}}(r,T)=-\frac{T}{12\pi}\sum_{j}\left[M_{j}^{3}(r)-m_{j}^{3}(r)\right]$
(2.18)

where $M_{j}$ are given by eigenvalues of the matrices
$\mathcal{M}_{j}+\Pi_{j}$ where $\mathcal{M}_{j}$ are mass matrices of bosons
and fermions and $\Pi_{j}$ are thermal self energies [51]. The summation of
$j$ includes all the bosons and fermions of the theory. Finally,

$V_{eff}^{total}(r,l,T)=V_{eff}(r,T)+V_{eff}(l,T)$ (2.19)

###  2.3 Two kinds of phase transitions

The finite temperature effective potential $V_{eff}$ in Eq. (2.19)
demonstrates a phase transition from the parity symmetric $v_{R}=v_{L}=0$
state at high temperature to the parity violating $v_{R}\neq v_{L}=0$ state at
low temperature. The transition occurs at a critical temperature $(T=T_{c})$
and can be of two types: first-order or second-order, depending on whether or
not the potential $V_{eff}$ has two degenerate global minima at $T_{c}$.
However for the left-right case, this classification needs to be modified for
two different reasons: finite causal horizon in the early Universe [68][69],
and the fact that there are two fields $\Delta_{L}$ and $\Delta_{R}$ obeying a
potential symmetric under $\Delta_{L}\leftrightarrow\Delta_{R}$. We may thus
identify two possibilities :

  * •

Causal horizon limited SOPT : SOPT is characterised by fluctuation length
$\xi=1/m_{r}=(\sqrt{2\rho_{1}}\eta)^{-1}$, the inverse mass of the Higgs
degree of freedom $r$ in Eq. (2.13). However, the phase transition can be
considered to have ended only after the temperature drops below the Ginzburg
temperature corresponding to the correlation length $\xi_{G}$

$\xi_{G}\simeq\frac{1}{2\rho_{1}\eta}$ (2.20)

This is determined by the condition that

$\xi_{G}^{3}\Delta_{c}V_{eff}=T_{G}$ (2.21)

which defines $T_{G}$ implicitly, where $\Delta_{c}V_{eff}$ is the difference
between free energy density in the two phases. In the early Universe, this
length is further restricted by the condition that the correlation length
cannot be diverging faster than the speed of light, which then implies that
the relevant correlation length is

$\xi_{causal}=\left(\frac{M_{Pl}}{\sqrt{\cal{N}}m_{r}^{2}T_{c}^{2}}\right)^{1/3}$
(2.22)

These considerations imply that as the correlation length reaches its maximum
value, uncorrelated regions with either $v_{L}\neq 0$ or $v_{R}\neq 0$ are
allowed to develop, with domain walls forming between them. This is unusual
for a SOPT but is true in the present theory. Thus $\xi_{causal}$ can be
considered to provide the scale of domain sizes just below the temperature
$T_{G}$.

  * •

Degenerate fields FOPT: In this case, while the true vacuum bubbles form
spontaneously, they can contain either $v_{L}\neq 0$ or $v_{R}\neq 0$. In
addition, as sketched in Fig. 1, after the merger of the same type of vacua,
there will arise a residual network of frustrated domain walls at a larger
length scale.

Figure 1: Cartoon explaining the emergence of residual domain walls after FOPT
in the case of two degenerate fields such as in the left-right symmetric
model. The symmetric false vacuum can be destabilized by the spontaneous
formation of bubbles of either L type vacuum with $v_{L}=0$ and $v_{R}\neq 0$
or R type with $v_{R}=0$ and $v_{L}\neq 0$. The domain wall sketched with
dotted lines eventually emerges as the bubbles of the same vacua merge to form
large homogeneous domains.

Using the tunneling probability formula in terms of the bounce solution, one
introduces a time scale: inverse duration of a phase transition ( see Eq.
(5.8)) $\beta$. So long as $\beta\gtrsim H,$ a sufficient number of bubbles
will be formed so that they can percolate[70, 35]. However, the percolation
can be either among the left-like or the right-like phases. Where the two
percolated regions meet, there is a domain wall, until the whole Universe is
filled with a frustrated network of domain walls. In the absence of any other
scales, it is reasonable to assume that over scales greater than the causal
horizon, the next near neighbour domain can randomly be in the same or the
other vacuum. This can lead to a frustrated network of scale larger than the
causal horizon by a factor $O(1)$ though not much large. For a space filling
network of irregular polyhedrons with an average of $S$ faces for each of the
polyhedrons, we may assume average enhancement by a factor$\sqrt{S}$ due to
merger with an average of $S/2$ near neighbours. The particle physics model,
therefore, needs to contain a mechanism for a slow disintegration of these
residual walls to not conflict with late cosmology.

The subsequent evolution of these residual domain walls depends on the bias
between the two types of approximately degenerate vacua. Finally, the false
domains will shrink out of existence according to the bias, producing GWs from
the late time evolution of the decay of these domain walls. Thus in addition
to the spontaneous bubbles of the usual FOPT, we get a contribution from the
CHL-SOPT type evolution of disintegrating DW. Since the two epochs are well
separated in time we can calculate each contribution individually.

We shall discuss the causal horizon limited SOPT in Section 3 and degenerate
fields FOPT in section 5. As we shall see, the contributions from the pure
FOPT and the CHL-SOPT are well separated in frequency space, allowing us to
piece together the result for the Degenerate Field FOPT case from the two
results.

##  3 The Causal Horizon Limited Second Order Phase Transition

###  3.1 GW from Domain Wall Decay

For the potential in Eq. (2.13), with two degenerate minima at $r=\pm\eta$,
one can find a static solution of the equation of motion with boundary
conditions that the two minima are realized at the extremities, i.e.
$r(x=\pm\infty)=\pm\eta$, of the form,

$r(x)=\eta\tanh\left(\frac{1}{2}\frac{x}{\delta_{\text{wall}}}\right)$ (3.1)

where $x$ is the spatial coordinate. This represents a domain wall extended
along the $x=0$ plane. $\delta_{\text{wall}}$ is the DW width, approximately
equal to the inverse of the mass of the field $r$,

$\delta_{\text{wall}}\sim m_{r}^{-1}=\left(\sqrt{2\rho_{1}}\eta\right)^{-1}$
(3.2)

Another important parameter of a domain wall is the tension, given by,

$\sigma_{\text{wall}}=\int_{-\infty}^{\infty}dx\rho_{r}=\frac{2}{3}m_{r}\eta^{2}$
(3.3)

where $\rho_{r}=\frac{1}{2}\left|\nabla r\right|^{2}+V_{0}(r)$ is the static
energy density of the field. These two parameters characterize a domain wall.

The evolution of the domain wall is determined by the Klein-Gordon equation.
As discussed in section 2.3, in case of a second order phase transition,
_frustrated_ domain walls are produced which move under the vacuum pressure
difference between the true and the false vacua. Following [53], we introduce
a small bias in the tree potential in Eq. (2.13) which breaks the degeneracy
of the two minima of the form,

$V_{\text{bias}}=\epsilon\eta r\left(\frac{r^{2}}{3}-\eta^{2}\right)$ (3.4)

where $\epsilon$ is a very small number. This term introduces an energy
difference between the two minima $\delta V=4\epsilon\eta^{4}/3$, which in
turn exerts a vacuum pressure on the domain wall $p_{V}\sim\delta V$. The
pressure difference should be large enough to make the walls disappear before
DWs start dominating the energy density of the Universe. The decay time is
given by,

$t_{\text{dec}}=\frac{1}{2}\sqrt{\frac{\rho_{1}}{2}}(\epsilon\eta)^{-1}$ (3.5)

However, due to the requirement that both vacua percolate, the energy
difference can’t be arbitrarily large. The ratio of the percolation
probabilities is given by $p_{+}/p_{-}\simeq e^{-4\delta V/\rho_{1}\eta^{4}}$
[34], which gives an upper bound on the bias
$\epsilon\leq-\frac{3}{16}\rho_{1}\log[(p_{+}/p_{-})_{\text{min}}]$. In a 3D
cubic lattice, if the probability of attaining a vacuum $p_{\pm}$ is greater
than a critical value $p_{c}=0.311$, then infinite clusters of the vacuum can
form [71]. Assuming $p_{+}=0.311$ and $p_{-}=1-0.311=0.689$ we get,

$\epsilon\leq 0.15\rho_{1}$ (3.6)

Figure 2: GW spectrum as seen today, from crumbling DWs for BP1, BP2, BP3 with
$v_{R}=10^{4},10^{5},10^{6}$ GeV respectively and the bias, $\epsilon=2\times
10^{-26}$, in comparison to _power-law-integrated-sensitivities_ of different
experiments. Figure 3: GW spectrum for BP2 with $\eta=10^{5}$ GeV and bias,
$\epsilon=10^{-27},10^{-26},10^{-25},10^{-24}$. The smaller the bias, the more
the amplitude but we can’t have an arbitrarily small bias since that will make
the domain walls enter the scaling regime. Figure 4: GW spectrum (in black)
for $(\eta,\epsilon)=(10^{10}\text{GeV},10^{-17})$ (solid),
$(\eta,\epsilon)=(10^{12}\text{GeV},10^{-13})$ (dashed) and
$(\eta,\epsilon)=(10^{13}\text{GeV},10^{-11})$ (dot-dashed). Here we have
taken $\rho_{1}=0.4$.

On the other hand, requirements for the well-known _scaling solution_ gives
another bound [72],

$\epsilon<4.7\times\left(\rho_{1}g_{*}\right)^{1/2}\frac{\eta}{M_{P}}$ (3.7)

where $g_{*}$ is the number of relativistic degrees of freedom at transition
temperature $T_{c}$ and $M_{P}$ is the Planck mass. Requiring that the walls
decay before they dominate the energy density, we get a lower bound on the
bias [72],

$\epsilon>\frac{16\pi}{9}\rho_{1}\left(\frac{\eta}{M_{P}}\right)^{2}$ (3.8)

###  3.2 Benchmark values

For evaluating the potential given in Eq. (2.14) and to see what type of phase
transitions take place, we use Mathematica. LRSM is very constrained in terms
of realization and parameters. In the triplet model, we set the $\beta_{i}$ s
to be 0 to obtain observable gauge and Higgs bosons [64]. We consider the
parameter space for our model while satisfying the necessary theoretical and
experimental constraints described in [51, 73]:

$\displaystyle v$
$\displaystyle=246\mathrm{GeV},v_{R}\in\left[10^{4},10^{6}\right]\mathrm{GeV},\tan\beta=\tan
10^{-3}$ $\displaystyle\lambda_{1}$
$\displaystyle=0.13,\lambda_{2}=0,\lambda_{3}\in[0,5],\lambda_{4}=0$
$\displaystyle\rho_{1}$
$\displaystyle\in[0,0.5],\rho_{2}\in[0,5],\rho_{3}\in[1,5],\rho_{4}=0$ (3.9)
$\displaystyle\alpha_{1}$
$\displaystyle=0,\alpha_{2}\in[0,0.5],\alpha_{3}\in[0,5]$
$\displaystyle\beta_{1}$ $\displaystyle=\beta_{2}=\beta_{3}=0$

We randomly generate sample points within these limits and check for the
global minima of the potential in Eq. 2.14 numerically using Mathematica’s
NMinimize function as well as using _CosmoTransitions_ [74] and match the two
results. Note that the temperature correction term in Eq. 2.2 as a function of
$r$, is not entirely smooth and may contain many local minima of negligible
depth; but they might confuse the NMinimize function while finding the global
minima of the effective potential. Therefore we approximate the effective
potential by a smooth polynomial function at different temperatures before
minimizing it. We look for both strong FOPT and SOPT in CosmoTransitions and
Mathematica. A continuous transition VEV from $r=0$ at high temperature to
$r=v_{R}(T)$ at low temperature would imply a second-order or a _cross-over_
phase transition whereas a barrier between two minima at the critical
temperature would imply a first-order phase transition. Out of all the
samples, we take 3 of the SOPT as benchmark points BP1, BP2, BP3 given in
Table 1 with energy scales $v_{R}=10^{4},10^{5},10^{6}$ GeV respectively for
the GW considerations from domain wall decay. In case of FOPT also, we choose
3 benchmark points named BPF1, BPF2, BPF3 given in Table 2 for calculating the
critical bubble profile in CosmoTransitions and using it to calculate the
relevant GW parameters in Mathematica as will be described in Sec. 5.

| BP1 | BP2 | BP3
---|---|---|---
$v/$GeV | 246 | 246 | 246
$v_{R}/$GeV | $10^{4}$ | $10^{5}$ | $10^{6}$
$\tan\beta$ | 0.001 | 0.001 | 0.001
$\lambda_{1}$ | 0.13 | 0.13 | 0.13
$\lambda_{2}$ | 0 | 0 | 0
$\lambda_{3}$ | 1.95024 | 0.66708 | 1.548816
$\lambda_{4}$ | 0 | 0 | 0
$\rho_{1}$ | 0.40055 | 0.411194 | 0.350746
$\rho_{2}$ | 0.88258 | 0.63584 | 0.847746
$\rho_{3}$ | 1.35815 | 1.711196 | 1.14436
$\rho_{4}$ | 0 | 0 | 0
$\alpha_{1}$ | 0 | 0 | 0
$\alpha_{2}$ | 0.27152 | 0.340437 | 0.2123735
$\alpha_{3}$ | 0.82883 | 0.603991 | 0.09274
$\beta_{1,2,3}$ | 0 | 0 | 0
$g$ | 0.65 | 0.65 | 0.65
$g_{\text{B-L}}$ | 0.4324 | 0.4324 | 0.4324
$y_{t}$ | 0.95 | 0.95 | 0.95
$y_{M}$ | 1 | 1 | 1
Table 1: Benchmark points for SOPT

Assuming the DWs vanish within the radiation-dominated era, the peak amplitude
and the peak frequency of the GWs produced, as seen at the present time, are
given by [75, 76],

$\displaystyle\Omega_{\mathrm{gw}}h^{2}\left(t_{0}\right)_{\text{peak }}$
$\displaystyle\simeq 5.20\times
10^{-20}\times\tilde{\epsilon}_{\mathrm{gw}}\mathcal{A}^{4}\left(\frac{10.75}{g_{*}}\right)^{1/3}\left(\frac{\sigma_{\text{wall
}}}{1\mathrm{TeV}^{3}}\right)^{4}\left(\frac{1\mathrm{MeV}^{4}}{\delta
V}\right)^{2}$ (3.10) $\displaystyle f_{\text{peak}}$ $\displaystyle\simeq
3.99\times
10^{-9}\mathcal{A}^{-1/2}\left(\frac{1\mathrm{TeV}^{3}}{\sigma_{\text{wall
}}}\right)^{1/2}\left(\frac{\delta
V}{1\mathrm{MeV}^{4}}\right)^{1/2}\mathrm{~{}Hz}$ (3.11)

where $t_{0}$ denotes the present time. $\mathcal{A}$ is the area parameter
$\simeq 0.8$ [76, 77] arising from $Z_{2}$ parity breaking of the LRSM, and
$\tilde{\epsilon}_{\text{gw}}$ is called the efficiency parameter $\simeq 0.7$
[76]. $g_{*}$ is the relativistic degrees of freedom at the time of the DW
decay, $g^{*}=134$ in our case. The spectrum is approximately given by,

$\left.\Omega_{\mathrm{GW}}\simeq\Omega_{\mathrm{GW}}\right|_{\text{peak
}}\times\left\\{\begin{array}[]{ll}\left(\frac{f_{\text{peak
}}}{f}\right)&\text{ for }f>f_{\text{peak }}\\\ \left(\frac{f}{f_{\text{peak
}}}\right)^{3}&\text{ for }f<f_{\text{peak }}\end{array}\right.$ (3.12)

We plot the GW spectrum for different benchmark points (as given in Table 1)
in Fig. 2 along with the power-law-integrated sensitivity curves (PLISC)s of
different GW experiments [78]. The 3 benchmark points (BP1, BP2, BP3) in Fig.
2 stands for $v_{R}=10^{4},10^{5},10^{6}$ GeV respectively with
$\epsilon=2\times 10^{-26}$. Instead of doing the traditional signal-to-noise
ratio (SNR) calculation of each signal against the noise profile of each
experiment to infer about the detectability of the signal [79, 80, 81], we
simply use the PLISCs to qualitatively state whether the signal will be
detectable or not in the future. A PLIS curve is calculated from the noise
profile of the experiment by considering the threshold SNR for the experiment.
The PLISCs given here are calculated assuming the threshold to be 1 and
observation time to be 1 year for interferometers and 20 years for PTA
experiments. Generally signal spectrum $\Omega_{\mathrm{signal}}$ is said to
be detectable in an experiment if the SNR for the signal is larger than the
threshold SNR for the experiment. But PLISCs provide a qualitative direct
graphical interpretation on the signal detectability without actually
calculating the SNR. A signal will be observable if
$\Omega_{\mathrm{signal}}(f)>\Omega_{\mathrm{PLIS}}(f)$ for some $f$. From
Fig. 2 it is clear that these signals are detectable at the PTA experiments
assuming a threshold SNR of 1. Fig. 3, shows the GW spectrum for BP2 with
different biases. Note that the parameter $\rho_{1}$ is primarily responsible
for determining the strength of the phase transition. Especially for bigger
values of $\rho_{1}\gtrsim 0.3$, the phase transition seems to be second order
or weak first order in nature, as can be seen in our benchmark points.
Considering this fact, we considered another set of benchmark points with
$\rho_{1}=0.4$ and
$(\eta,\epsilon)=(10^{10}\text{GeV},10^{-17}),(10^{12}\text{GeV},10^{-13}),(10^{13}\text{GeV},10^{-11})$
which should give detectable GW peaks at the interferometers as seen in Fig.
4. Note that we did not use the peak-integrated sensitivity curves (PISC)s
here which also encode the SNR information of the signal, because PISCs are
not available for the PTA experiments [78].

##  4 Numerical Simulation

###  4.1 Equation of Motion

We consider the tree level potential as in Eq. (2.13) to simulate domain wall
decay (see [53]). It is symmetric under $r\rightarrow-r$. It has two
degenerate minima, at $\pm\eta$. The height of the potential barrier is
$\Delta V\approx\rho_{1}\eta^{4}/4$. The evolution of the field in an
expanding FLRW metric is given by the Klein-Gordon equation,

$\ddot{r}+H\dot{r}-\frac{\nabla^{2}}{a^{2}}r+\frac{dV}{dr}=0$ (4.1)

where dot denotes derivative w.r.t. $t$ and $a=a(t)$ is the scale factor. This
equation also describes the evolution of domain walls.

We use the notation of [53] and in the appendix A we summarize the methodology
of the calculation, which utilizes the Green’s function method. The energy
density of GW is given by Eq. (A.2). Replacing the ensemble average by an
average over the volume of the comoving box, we obtain,

$\rho_{\mathrm{gw}}=\frac{1}{32\pi
Ga^{4}}\frac{1}{V}\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}\bar{h}_{ij}^{\prime}(\tau,\mathbf{k})\bar{h}_{ij}^{\prime*}(\tau,\mathbf{k})$
(4.2)

Substituting the form of $\bar{h}_{ij}$ from Eq. (A.12) and (A.2), upto first
order of $aH$, we obtain,

$\displaystyle\rho_{\mathrm{gw}}=\frac{2\pi^{2}G}{a^{4}V}\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}\frac{1}{k^{2}}\times\sum_{ij}$
$\displaystyle\left\\{\left|\int_{x_{i}}^{x_{f}}dx^{\prime}\sqrt{x^{\prime}}a\left(x^{\prime}\right)N_{\nu}\left(x^{\prime}\right)T_{ij}^{\mathrm{TT}}\left(\tau^{\prime},\mathbf{k}\right)\right|^{2}\right.$
$\displaystyle\left.+\left|\int_{x_{i}}^{x_{f}}dx^{\prime}\sqrt{x^{\prime}}a\left(x^{\prime}\right)J_{\nu}\left(x^{\prime}\right)T_{ij}^{\mathrm{TT}}\left(\tau^{\prime},\mathbf{k}\right)\right|^{2}\right\\}$
(4.3)

where we have used approximations for $k\tau\gg 1$,

$J_{\nu}(k\tau)\rightarrow\sqrt{\frac{2}{\pi(k\tau)}}\cos\left(k\tau-\frac{\nu\pi}{2}-\frac{\pi}{4}\right),\quad
N_{\nu}(k\tau)\rightarrow\sqrt{\frac{2}{\pi(k\tau)}}\sin\left(k\tau-\frac{\nu\pi}{2}-\frac{\pi}{4}\right)$
(4.4)

and averaged over a period of oscillation of $\bar{h}_{ij}(\tau,\mathbf{k})$
with time. The fraction of energy density of GWs given by Eq. (A.15) becomes,

$\displaystyle\Omega_{\mathrm{gw}}(k,t)=\frac{2G^{2}k}{3Va(t)^{4}H(t)^{2}}\int
d\Omega_{k}$
$\displaystyle\times\sum_{ij}\left\\{\left|\int_{x_{i}}^{x_{f}}dx^{\prime}\sqrt{x^{\prime}}a\left(x^{\prime}\right)N_{\nu}\left(x^{\prime}\right)T_{ij}^{\mathrm{TT}}\left(\tau^{\prime},\mathbf{k}\right)\right|^{2}\right.$
$\displaystyle\left.\quad+\left|\int_{x_{i}}^{x_{f}}dx^{\prime}\sqrt{x^{\prime}}a\left(x^{\prime}\right)J_{\nu}\left(x^{\prime}\right)T_{ij}^{\mathrm{TT}}\left(\tau^{\prime},\mathbf{k}\right)\right|^{2}\right\\}$
(4.5)

where $\Omega_{k}$ is a unit vector representing the direction of $\mathbf{k}$
and $d\Omega_{k}=d\cos\theta dr$. The TT part of the stress-energy tensor is
computed by applying the projection operator in the momentum space,

$\displaystyle T_{ij}^{\mathrm{TT}}(\tau,\mathbf{k})$
$\displaystyle=\Lambda_{ij,kl}(\hat{k})T_{ij}(\tau,\mathbf{k})$
$\displaystyle=\Lambda_{ij,kl}(\hat{k})\left\\{\partial_{k}r\partial_{l}r\right\\}(\tau,\mathbf{k})$
(4.6) $\displaystyle\Lambda_{ij,kl}(\hat{k})$
$\displaystyle=P_{ik}(\hat{k})P_{jl}(\hat{k})-\frac{1}{2}P_{ij}(\hat{k})P_{kl}(\hat{k})$
(4.7) $\displaystyle P_{ij}(\hat{k})$
$\displaystyle=\delta_{ij}-\hat{k}_{i}\hat{k}_{j}$ (4.8)

where $\hat{k}=\mathbf{k}/|\mathbf{k}|$, and
$\left\\{\partial_{k}r\partial_{l}r\right\\}(\tau,\mathbf{k})$ is the Fourier
transform of $\partial_{k}r(\tau,\mathbf{x})\partial_{l}r(\tau,\mathbf{x})$.

###  4.2 Simulation Setup

The simulation takes place in a cuboid of size $L_{x}=L_{y}=L_{z}=b$ in terms
of $\eta^{-1}$, discretized into $N^{3}=256^{3}$ lattice points. Our
simulation box is a comoving box, which also expands with the expansion of the
universe. The problem with this setup is that the dynamic range of the
simulation is limited. This is because we want our comoving box to be bigger
than the correlation length, and also want to resolve the domain wall width,
at all times during the simulation. Domain wall width remains constant in real
space, i.e. decreases as $a(t)^{-1}$ within the simulation box.

Let our physical lattice spacing be $\delta x_{phy}=a(t)\times b/N=a(t)\times
50/256$. Assuming initial time of simulation $t_{i}=1$ and setting
$a(t_{i})=1$, the ratio of domain wall width to lattice spacing is,

$\frac{\delta_{w}}{\delta
x_{phy}}=\frac{N}{b\rho_{1}^{1/2}}\left(\frac{t}{t_{i}}\right)^{-\gamma}$
(4.9)

where $\gamma$ is from Eq. (A.6). In our case, the entire process occurs in
the radiation dominated universe ($\gamma=1/2$), and we choose $\rho_{1}=0.4$,
$b=50$ and vary $\epsilon$ from $0$ to $0.015$ for the simulations.
Corresponding to each $\epsilon$, we get an approximate wall decay time
$t_{\text{dec}}$ from Eq. (3.5). But we see that for $t>61$ (in units of
$\eta^{-1}$) the above ratio is $>1$, which means that the walls are not
resolvable anymore. This limits the upper bound on the time integral of the GW
calculation in Eq. (4.1) to be the resolution time $t=t_{\text{res}}=61$, i.e.
we consider the source to be active from $t=1$ to $t=61$. Beyond that, we will
get unphysical results. We chose the maximum time of simulation to be $151$ to
see how the spectrum changes beyond $t=61$ and plot them in Figures 11 and 12.
The different plots in each graph correspond to what upper limit we take in GW
calculation in Eq. (4.1). Ideally $t_{\text{dec}}$ should be smaller than 61
in our simulation to obtain realistic results but we also consider larger
$t_{\text{dec}}$ and interpret the results carefully as discussed in Sec. 4.4.

The correlation length determines the initial size of the domains. Assuming
$m_{r}$ and $T_{c}$ to be of the order of $\eta$, dimensionally we can find
$\xi_{\text{causal}}/\eta^{-1}\sim\left(\frac{M_{\text{pl}}}{\eta}\right)^{1/3}$,
which is around 10 for $\eta=10^{16}$ GeV and 1000 for $\eta=10^{10}$ GeV. So
we take $\eta=10^{16}$ GeV for our simulation so that the domain sizes are of
the order $\mathcal{O}(10)$.

Spectrum calculation procedure:

  1. 1.

Simulate the field for a specific time step to obtain $r(\tau,\mathbf{x})$.

  2. 2.

Calculate $\left\\{\partial_{k}r\partial_{l}r\right\\}(\tau,\mathbf{k})$ by
taking Discrete Fourier Transform (DFT) of
$\partial_{k}r(\tau,\mathbf{x})\partial_{l}r(\tau,\mathbf{x})$

  3. 3.

Calculate $T_{ij}^{\mathrm{TT}}(\tau,\mathbf{k})$ using Eq. (4.1)-(4.8).

  4. 4.

Perform the time integration in Eq. (4.1) and the sum over $i,j$. We used a
simple Simpson’s 3/8th rule for the time integration.

  5. 5.

Perform the angular integration $\int d\Omega_{k}$ using Monte Carlo method.
Multiply by the appropriate constant as in Eq. (4.1).

####  4.2.1 Initial Field Configuration

We assume that domain walls are already produced at $t_{i}$ with a
characteristic length scale given by the correlation length, motivated by Eq.
(2.22). So we take a Sin function with the wavelength corresponding to twice
the correlation length at the initial time as our dominant mode, i.e.
$k_{\text{dominant}}=\pi/\xi_{\text{causal}}$. One wavelength of the Sin
function consists of two domains of opposite vacua. We randomize the
configuration by adding 4-5 more modes to it with small coefficients. Finally,
we take a Tanh of this series of Sin functions which squeezes the rise and
fall of the Sin curves, as a good model of the domain wall kink. The final
results are obtained as broad dependence of the power on frequency and are
insensitive to minor changes in the initial conditions. The needed
configuration corresponds to a mesh of cubical domains where each domain is
surrounded by the opposite kind of domains and can be created by defining the
initial field values as,

$r(x,y,z)=\tanh{\left[N_{\text{norm}}\sum_{n=1}^{i=5}C_{n}\sin{\left(k_{n}x+\phi_{i}\right)}\sin{\left(k_{n}y+\psi_{i}\right)}\sin{(k_{n}z+\chi_{i})}/\delta_{\text{wall}}\right]}$
(4.10)

where $N_{\text{norm}}$ is a normalizing factor, $C_{1}=1$ and the rest of the
$C_{i}$ are small numbers of the order $\mathcal{O}(10^{-1})$. $k_{n}=\pi
n/\xi_{\text{causal}}$ represents the initial characteristic comoving
wavenumber. Every dimensional parameter is scaled by $\eta$. Finally
$\phi_{i},\psi_{i},\chi_{i}$ are random phases. It turns out that the
resulting spectrum depends mostly on the dominant fluctuation scales rather
than the exact initial configuration.

In Figure 9, the evolution of the domain wall network is shown for
$\epsilon=0$ and $\epsilon=0.007$. The effect of the bias is seen as false
domains reduce in size and domain walls decay, while there’s no preferred
domain in the zero bias case.

(a) $t=21,\epsilon=0$

(b) $t=21,\epsilon=0.007$

(a) $t=51,\epsilon=0$

(b) $t=51,\epsilon=0.007$

Figure 6:

(a) $t=81,\epsilon=0$

(b) $t=81,\epsilon=0.007$

Figure 7:

(a) $t=111,\epsilon=0$

(b) $t=111,\epsilon=0.007$

Figure 8:  (a) $t=141,\epsilon=0$ Figure 9: Evolution of domain walls for
$\epsilon=0$ on the left column and $\epsilon=0.007$ on the right column.
Image (i) is the final form of the no-bias case after the simulation.

###  4.3 Scaling Regime

An important feature of domain wall evolution without bias is the era of
scaling when the characteristic length scales like the radius of curvature of
domain walls or domain size is comparable to the Hubble radius,

$R\sim H^{-1}\sim t$ (4.11)

This phase of evolution is found both numerically [82, 83, 84, 85, 86, 87] and
analytically [88, 89, 90] in the literature, and it corresponds to “at least
one domain wall per horizon" scenario. While in the scaling regime, the energy
density of the domain walls is given by,

$\rho_{\text{wall}}\sim\sigma_{\text{wall}}R^{2}/R^{3}\sim\sigma_{\text{wall}}/t$
(4.12)

with $\sigma_{\text{wall}}$ given in Eq. (3.3), Eq. (4.12) can be expressed
as,

$A/V\propto\tau^{-1}$ (4.13)

where $A/V$ is the comoving area density of the walls and $\tau$ is the
conformal time. To calculate this quantity we define a wall to be the surface
where the field value is 0. Therefore, to measure the surface area of walls in
the 3D lattice, first we have to figure out where the field crosses zero value
within the lattice. To do this, let us call two lattice points that are
nearest neighbours as links. These links are oriented in $x,y,z$ directions
and we measure the signs of the field at the two ends of each link. If the
signs are different at the two ends, we have found a zero crossing of the
field. For all such possible links at all the lattice points within the cube,
we do the same and assign a variable $\delta_{\pm}$ to each link, which takes
on value $+1$ if we found a zero crossing, otherwise $0$. Once we find such a
link that penetrates a wall, we add one unit grid area (area of the square
made by nearest points) to the total surface area of walls but divided by a
weight factor. The weight factor must be the average number density of points
on the surface, based on the orientation of the surface at the location of the
link. The orientation can be found through the unit normal vector to the wall
(given by the gradient of the field) at that location. Finally, we sum all
such contributions from all the links with $\delta_{\pm}=1$ to obtain the
total wall area [82],

$A=\Delta A\sum_{\text{all
links}}\delta_{\pm}\frac{|\nabla\phi|}{|\partial_{x}\phi|+|\partial_{y}\phi|+|\partial_{z}\phi|}$
(4.14)

where $\Delta A$ is the unit grid area. In our case $\Delta A=\Delta
x^{2}=\Delta y^{2}=\Delta z^{2}=(b/N)^{2}$. We calculate this quantity divided
by volume $b^{3}$ of the cube at different times and plot in Fig. 10. The
walls seem to enter the scaling regime after $\tau\sim 10$.

(a) Area density with conformal time.

(b) Mean field value with time.

Figure 10: Evolution of a) comoving area density with conformal time and b)
mean field value with $t$. In a), for the no-bias case the walls seem to
follow scaling solutions after $\tau\sim 10$ where area density varies as
$\tau^{-1}$ as mentioned in Sec. 4.3. In b) we verify the equation for decay
time (3.5) upto $\mathcal{O}(1)$ by measuring the mean field value with time
for different simulations. The mean crosses the $r=1/\sqrt{2}$ line at around
$t=t_{\text{dec}}$.

###  4.4 Gravitational Wave Spectrum

We normalize the amplitude of GW given in Eq. (4.1) by,

$\Omega_{\eta}\equiv\frac{\rho^{\eta}_{gw}}{\rho(t_{i})}=\frac{8\pi}{3\beta^{2}}G^{2}\eta^{4}$
(4.15)

where $\rho^{\eta}_{gw}=G\eta^{6}$ is the energy density of a source of
characteristic scale $\eta$ estimated from Quadrupole approximation, and
$\rho(t_{i})=3H^{2}(t_{i})/8\pi G=3\beta^{2}\eta^{2}/8\pi G$ is the critical
energy density at the initial time. We plot the energy spectrum for the no-
bias case in Fig. 11 and with different biases in Fig. 12. The horizontal axis
is the amplitude of the comoving wave vector. It is related to the frequency
as in Eq. (A.16). Notice that in the Monte Carlo method used to perform the
angular integration $d\Omega_{k}$, we define the shell of radius $k$ and
average out the values of the fields at the lattice points within that shell.
Therefore in the small $k$ region, the number of points within the shell is
also small which generates a statistical error in the calculation, i.e. the
lower $k$ range of $\mathcal{O}(0.1)$, where the peak occurs, is not properly
resolvable. We plot the error bars for the Monte Carlo integration to show the
statistical uncertainties in the calculation. Other sources of error include
numerical errors from the usage of a simple difference method for simulation
and Simpson’s rule for time integration. Therefore it is difficult to provide
conclusive remarks about the spectrum with our small dynamical range. However,
we can deduce some information about the spectrum as follows:

  1. 1.

As discussed in Sec. 4.2 we fit the curve (solid black line) at $t\sim
t_{\text{res}}\sim 61$ for different cases. We understood that beyond
$t_{\text{res}}$, we will get unphysical results. Those unphysical results
manifest in the graphs as overall higher amplitude if
$t_{\text{dec}}>t_{\text{res}}$, and resonance, seen as a bump around
$k\sim\mathcal{O}(10)$, corresponding to the wall width. The reason for the
resonance is that the wall width becomes equivalent or greater than the
physical lattice spacing beyond $t_{\text{res}}$. Therefore we will consider
the fitted line to be the actual spectrum from the simulation.

  2. 2.

It is seen that the spectrum peaks at around $k\sim 0.4-0.5$ which corresponds
to the Hubble size at the time of GW production, $k_{h}/a(t)=2\pi H(t)$ at
$t\sim 61$, as seen from the fitted line.

  3. 3.

The spectrum shows a dependence of $k^{-1}$ beyond the peak as seen in Fig. 11
for the no-bias case. If we evolve the field with a bias explicitly, an
enhancement towards the high frequency region is seen which can be interpreted
as contributions coming from small-sized walls as the simulation proceeds. In
the no-bias case, walls don’t crumble to give such a contribution. For
observation in experiments, such high frequency enhancement is irrelevant
since we are interested only in the peak at $k_{h}$. Thus we conclude that the
spectrum varies as $k^{-1}$ for $k>k_{h}$.

  4. 4.

Now due to the high error in the small $k$ range, the behavior of the spectrum
in that region is unclear. Therefore we assume a $k^{3}$ dependence for
$k<k_{h}$ considering requirements for causality [91]. These conclusions match
with [75, 76], therefore in Sec. 3, we use the same general formulas as given
in them for the GW spectrum generated at arbitrary scale with arbitrary bias,
as seen today. The formulas are given in (3.10) and (3.11).

Figure 11: Normalized GW spectrum vs comoving wavenumber with no bias.
Assuming the walls suddenly vanish at different $t$, we can calculate the bias
for the corresponding graph from Eq. (3.5) as shown in the legend. Error is
calculated as $\sigma/\sqrt{N}$ and shown by error bars. The spectrum peaks
around $k\sim k_{h}\sim\mathcal{O}(0.4-0.5)$, corresponding to the Hubble
radius for $t\sim 61$. The black dashed line is a curve fit at
$t=t_{\text{res}}=61$ which shows a dependence of $k^{-1}$ at the high
frequency region. Note that high errors in the low frequency region make it
difficult to exactly see the dependence of the spectrum for small $k$.
Therefore from causality requirements [91] we assume a $k^{3}$ dependence at
$k<k_{h}$. Beyond $t_{\text{res}}$, we get unphysical results as discussed in
Sec. 4.2.

(a) $\epsilon=0.0015$

(b) $\epsilon=0.003$

(c) $\epsilon=0.007$

(d) $\epsilon=0.015$

Figure 12: Spectrum for different biases. The wall decay time $t_{\text{dec}}$
for each case is calculated using the formula (3.5). The curve fit line has a
smaller slope compared to that in Fig. 11, which shows an enhanced
contribution in the high frequency region as compared to the no-bias case. The
more the bias more the enhancement and at $\epsilon=0.015$, the spectrum is
almost flat between the Hubble frequency and wall width frequency around
$k\sim\mathcal{O}(10)$. This enhancement comes from the small-sized walls
towards the end of the simulation.

##  5 The Degenerate Fields First Order Phase Transition

In this case, the phase transition begins with the usual spontaneous bubble
formation and mutual collisions but towards the end, leaves behind two-field
domain walls of the type discussed in [42, 43]. These should disintegrate on
their own, giving an additional contribution to the GW. We begin with the
contributions from a standard FOPT for our model, since the additional
contribution due to the residual DW of the Degenerate Field FOPT arises
independently. Gravitational wave spectrum from FOPT is characterized by two
parameters $\alpha$ and $\beta$, corresponding to the strength of the phase
transition and the inverse rate of tunneling. In this section, we will give
relevant definitions of these two parameters and speculate the observability
of the GW spectrum for typical values for them.

First order phase transition (FOPT) takes place through the nucleation of
bubbles within which $r$ takes the true VEV, $v_{R}(T)$. The field $r$ tunnels
from a metastable state $v_{R}(T)=0$ to the stable ground state $v_{R}(T)\neq
0$ considering enough supercooling. These bubbles are spherically symmetric
bounce solutions [92]. After nucleation, these bubbles of true VEV expand
through the plasma in the expanding metastable Universe, and depending on the
strength of interaction of the bubble wall with the plasma, either the bubble
walls accelerate to attain a terminal velocity before collisions take place
(non-runaway scenario in plasma) [93], or they keep accelerating without bound
to reach approximately the speed of light (runaway scenario) [93, 94, 95, 96].
There is a third scenario for super-strong phase transitions, where plasma
effects can be ignored and they attain speed of light (runaway in vacuum)
[93]. The nucleation rate of bubbles at temperature $T$ is given by [51, 44,
45, 97],

$\Gamma(T)\simeq T^{4}\left(\frac{S_{3}(T)}{2\pi
T}\right)^{\frac{3}{2}}\mathrm{e}^{-S_{3}(T)/T}$ (5.1)

where $S_{3}(T)$ is the three-dimensional Euclidean action evaluated for the
bounce solution which is obtained by solving the equation of motion,

$\frac{\mathrm{d}^{2}r}{\mathrm{~{}d}x^{2}}+\frac{2}{x}\frac{\mathrm{d}r}{\mathrm{~{}d}x}=\frac{\mathrm{d}V_{\mathrm{eff}}(r,T)}{\mathrm{d}r}$
(5.2)

with boundary conditions $dr/dx=0$ at $x=0$ and $r\rightarrow 0$ as
$x\rightarrow\infty$ where $x$ denotes the 3D radial coordinate. We use
CosmoTransitions to solve this equation and to get the critical bubble profile
$r(x)$. Once the profile is known we can calculate $S_{3}(T)$ by,

$S_{3}(T)=4\pi\int_{0}^{\infty}dxx^{2}\left[\frac{1}{2}\left(\frac{dr}{dx}\right)^{2}+V_{\text{eff}}(r(x),T)\right]$
(5.3)

Next, the nucleation temperature $T_{n}$ is the temperature at which there’s
at least one bubble created per horizon volume,

$\int_{T_{n}}^{T_{c}}\frac{\mathrm{d}T}{T}\frac{\Gamma(T)}{H(T)^{4}}=1$ (5.4)

where the Hubble parameter is given by,

$H(T)^{2}=\frac{\rho_{\mathrm{rad}}(T)+\rho_{\mathrm{vac}}(T)}{3\mathcal{M}_{\mathrm{Pl}}^{2}}=\frac{1}{3\mathcal{M}_{\mathrm{Pl}}^{2}}\left(\frac{\pi^{2}}{30}g_{*}T^{4}+\Delta
V(T)\right)$ (5.5)

with $g_{*}=134$ is the relativistic degrees of freedom for our model and
$\mathcal{M}_{pl}=2.435\times 10^{18}$ GeV is the reduced Planck Mass. Vacuum
energy density $\rho_{\mathrm{vac}}(T)$ is calculated from $\Delta
V(T):=V_{\mathrm{eff}}(0,T)-V_{\mathrm{eff}}\left(v_{R}(T),T\right)$. The
condition (5.4) can be approximated by $\Gamma(T)\sim H(T)^{4}$ and assuming
$\Gamma(T)\sim T^{4}e^{-S_{3}(T)/T}$ and radiation domination, i.e $H(T)\sim
1.66\sqrt{g_{*}}T^{2}/3M_{\text{Pl}}$ with $M_{\text{Pl}}=1.22\times 10^{19}$
GeV, we get,

$\frac{S_{3}(T_{n})}{T_{n}}\sim-4\log\frac{1.66\sqrt{g_{*}}T_{n}}{M_{\text{Pl}}}\sim
100-130$ (5.6)

for $g_{*}=134$ and $T_{n}\sim 10^{6}-10^{4}\text{ GeV}$. Once we have
$T_{n}$, we can define the parameters [49, 98],

$\alpha=\frac{1}{\rho_{\mathrm{rad}}\left(T_{n}\right)}\left(\Delta
V\left(T_{n}\right)-\left.\frac{T_{n}}{4}\frac{\partial\Delta V(T)}{\partial
T}\right|_{T=T_{n}}\right)$ (5.7)

which characterizes the strength of the FOPT; and the inverse duration of the
transition given by,

$\beta=\left.H\left(T_{n}\right)T_{n}\cdot\frac{\mathrm{d}\left(S_{3}/T\right)}{\mathrm{d}T}\right|_{T=T_{n}}$
(5.8)

For strong FOPT, the typical values for these two parameters for our model
are, $\alpha\in[0.001,0.1]$ and $\beta/H_{*}\in[10^{2},10^{4}]$ [73] where
$H_{*}=H(T_{n})$. Alternatively the condition $v_{R}(T_{c})/T_{c}>1$ where
$v_{R}(T)$ is the temperature dependent VEV of the field $r$ is also used to
characterize a strong FOPT. As mentioned in Sec. 3.2 this condition is
satisfied primarily for a small value of $\rho_{1}$ [51]. Therefore, in Table
2, we take 3 benchmark points with small $\rho_{1}$ values which signals
strong FOPT, which we will use to plot the GW spectrum. Apart from the
critical and nucleation temperature it is necessary to estimate the percoation
temperature $T_{p}$, defined as the temperature at which a significant portion
of the Hubble volume is attains true vacuum and percolation takes place, and
make sure that it is not far from the nucleation temperature. Following [50],
we calculate the percolation temperature by defining the probability of a
point still in false vacuum at temperature $T$ to be $P(T)=e^{-I(T)}$ where
$I(T)$ is the volume of true vacuum per unit comoving volume given by [99],

$I(T)=\frac{4\pi}{3}\int_{T}^{T_{c}}dT^{\prime}\frac{\Gamma(T^{\prime})}{T^{\prime
4}H(T^{\prime})}\left(\int_{T}^{T^{\prime}}\frac{d\tilde{T}}{H(\tilde{T})}\right)^{3}$
(5.9)

which under the assumption
$H(T)\sim\frac{1}{\sqrt{3}\mathcal{M}_{\text{pl}}}\frac{\pi^{2}\sqrt{g_{*}}}{\sqrt{30}}T^{2}$
takes the form,

$I(T)=\frac{4\pi}{3T^{3}}\left(\sqrt{3}\mathcal{M}_{\text{pl}}\frac{\sqrt{30}}{\pi\sqrt{g_{*}}}\right)^{4}\int_{T}^{T_{c}}dT^{\prime}\frac{\Gamma(T^{\prime})}{T^{\prime
6}}\left(1-\frac{T}{T^{\prime}}\right)^{3}$ (5.10)

The percolation temperature is calculated using $I(T_{p})=0.34$ [100]. As long
as $T_{p}$ is of the same order as $T_{n}$ we can assume that bubbles
nucleated at $T_{n}$ will percolate and take $T_{n}$ to be the proxy for
temperature in the calculation of $\alpha$ and $\beta$. The results of our
calculations for the 3 benchmark points BPF1, BPF2, BPF3 are shown in Table 3.

Now, there are three sources of GWs in case of a first-order phase transition;
i) bubble wall collision [101], ii) sound waves [91] and iii)
magnetohydrodynamic (MHD) turbulence [102]. The spherical symmetry of the
nucleated bubble is broken when the bubbles collide, which is the basic
requirement for GW production. Thus, the total GW spectrum is given by,

$h^{2}\Omega_{\mathrm{GW}}\simeq
h^{2}\Omega_{r}+h^{2}\Omega_{\mathrm{sw}}+h^{2}\Omega_{\mathrm{turb}}$ (5.11)

The first term is the contribution from bubble collisions which can be
calculated from the field evolution. The second and third terms are the sound
wave. In order to calculate the contributions from each of these three
sources, we need to distinguish between the three scenarios mentioned above,
namely nonrunaway, runaway and runaway in vacuum. To do this we compare our
$\alpha$ parameter with the $\alpha_{\infty}$ parameter [93] and found that
$\alpha$ is always smaller than $\alpha_{\infty}$, which corresponds to the
nonrunaway scenario. In this scenario the bubble walls expand for enough time
to attain terminal velocity before collisions. Therefore the contribution to
GW spectrum from collisions is negligible compared to the contribution from
sound waves and turbulences. This is because the energy released in collisions
comes from the gradient of the field $r$ which scales as the surface area of
the bubble, while the energy from sound waves and turbulences come from the
bulk of the plasma which scale as the volume.

###  5.1 Sound Waves

The movement of the DWs through the plasma creates pressure waves in the
plasma. The contribution of such sound waves to GW is given from numerical fit
[3, 103, 104],

$h^{2}\Omega_{\mathrm{sw}}(f)=2.65\times
10^{-6}\left(\frac{H_{*}}{\beta}\right)\left(\frac{\kappa_{sw}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{\frac{1}{3}}v_{w}\left(\frac{f}{f_{\mathrm{sw}}}\right)^{3}\left(\frac{7}{4+3\left(f/f_{\mathrm{sw}}\right)^{2}}\right)^{7/2}\Upsilon$
(5.12)

where $v_{w}$ is the wall velocity which we assume as [105],

$v_{w}=\begin{cases}\sqrt{\frac{\Delta
V(Tn)}{\alpha\rho_{\text{rad}}}},&\text{if }\sqrt{\frac{\Delta
V(Tn)}{\alpha\rho_{\text{rad}}}}\geq v_{J}\\\ 1,&\text{otherwise}\end{cases}$
(5.13)

where $v_{J}$ is the Jouguet velocity,
$v_{J}=\frac{1}{\sqrt{3}}\frac{1+\sqrt{3\alpha^{2}+2\alpha}}{1+\alpha}$. The
efficiency factor $\kappa_{sw}$ is given by [78, 93],

$\kappa_{\mathrm{sw}}=\begin{cases}\frac{\alpha}{0.73+0.083\sqrt{\alpha}+\alpha},&\text{if
}v_{w}\geq v_{w}^{\alpha}\\\ \frac{6.9\alpha
v_{w}^{6/5}}{1.36-0.037\sqrt{\alpha}+\alpha},&\text{if
}v_{w}<v_{w}^{\alpha}\end{cases}$ (5.14)

where,

$v_{w}^{\alpha}=\left[\frac{1.36-0.037\sqrt{\alpha}+\alpha}{6.9(0.73+0.083\sqrt{\alpha}+\alpha)}\right]^{5/6}$
(5.15)

Lastly the factor $\Upsilon=1-\frac{1}{\sqrt{1+2\tau_{sw}H_{*}}}$ is a
suppression factor dependent on lifetime of sound waves $\tau_{sw}$ [104]. It
can be parameterized by writing $\tau_{sw}\sim R_{*}/\overline{U}_{f}$, where
$R_{*}=(8\pi)^{1/3}v_{w}/\beta$ and
$\overline{U}_{f}=\sqrt{3\kappa_{\mathrm{sw}}\alpha/4}$ represents mean bubble
distance and root-mean-squared fluid velocity [106]. The peak frequency is
given by,

$f_{\mathrm{sw}}=\frac{1.9\times
10^{-5}}{v_{w}}\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100\mathrm{GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}}\mathrm{Hz}$
(5.16)

###  5.2 MHD Turbulence

Turbulent motion of the fully ionized plasma also contributes to the GW
spectrum, modeled by [3, 107, 108],

$h^{2}\Omega_{\mathrm{turb}}(f)=3.35\times
10^{-4}\left(\frac{H_{*}}{\beta}\right)\left(\frac{\kappa_{\text{turb
}}\alpha}{1+\alpha}\right)^{\frac{3}{2}}\left(\frac{100}{g_{*}}\right)^{1/3}v_{w}\frac{\left(f/f_{\text{turb
}}\right)^{3}}{\left[1+\left(f/f_{\text{turb
}}\right)\right]^{\frac{11}{3}}\left(1+8\pi f/h_{*}\right)}$ (5.17)

where $\kappa_{\mathrm{turb}}=0.05\kappa_{\mathrm{turb}}$ [93] and,

$h_{*}=16.5\cdot
10^{-6}\left(\frac{T_{n}}{100\mathrm{GeV}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}\mathrm{~{}Hz}$
(5.18)

The peak frequency is given by,

$f_{\text{turb }}=\frac{2.7\times
10^{-5}}{v_{w}}\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100\mathrm{GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}}\mathrm{Hz}$
(5.19)

| BPF1 | BPF2 | BPF3
---|---|---|---
$v/$GeV | 246 | 246 | 246
$v_{R}/$GeV | $10^{4}$ | $10^{5}$ | $10^{6}$
$\tan\beta$ | 0.001 | 0.001 | 0.001
$\lambda_{1}$ | 0.13 | 0.13 | 0.13
$\lambda_{2}$ | 0 | 0 | 0
$\lambda_{3}$ | 3.452 | 1.3374 | 2.284
$\lambda_{4}$ | 0 | 0 | 0
$\rho_{1}$ | 0.0264 | 0.0386 | 0.0787
$\rho_{2}$ | 1.0542 | 0.5973 | 0.2495
$\rho_{3}$ | 2.6854 | 2.0801 | 3.3698
$\rho_{4}$ | 0 | 0 | 0
$\alpha_{1}$ | 0 | 0 | 0
$\alpha_{2}$ | 0 | 0 | 0
$\alpha_{3}$ | 0.5796 | 0.4494 | 0.3537
$\beta_{1,2,3}$ | 0 | 0 | 0
$g$ | 0.65 | 0.65 | 0.65
$g_{\text{B-L}}$ | 0.4324 | 0.4324 | 0.4324
$y_{t}$ | 0.95 | 0.95 | 0.95
$y_{M}$ | 0.9362 | 0.7875 | 0.9255
Table 2: Benchmark points for FOPT | BPF1 | BPF2 | BPF3
---|---|---|---
$T_{c}/$GeV | 2913 | 28829 | 323004
$T_{n}/$GeV | 1650 | 21020 | 275040
$T_{p}/$GeV | 1557.2 | 20426 | 269050
$v_{R}(T_{c})/T_{c}$ | 4.37146 | 3.9021 | 3.04868
$\alpha$ | 0.395862 | 0.141741 | 0.025492
$\alpha_{\infty}$ | 0.789406 | 0.334233 | 0.0722176
$\beta/H_{*}$ | 200.641 | 547.352 | 648.779
$\sqrt{\frac{\Delta V(Tn)}{\alpha\rho_{\text{rad}}}}$ | 0.927674 | 0.82288 | 0.391911
$v_{J}$ | 0.878237 | 0.802155 | 0.692529
$v_{w}^{\alpha}$ | 0.275777 | 0.303166 | 0.325535
Table 3: Phase transition parameters calculated for the benchmark points for
FOPT given in Table 2. Here $\sqrt{\frac{\Delta
V(Tn)}{\alpha\rho_{\text{rad}}}}>v_{J}$ implies we take $v_{w}=1$ for BPF1 and
BPF2 while for BPF3 $v_{w}=0.391911$ as in Eq. (5.13).  Figure 13: GW spectrum
as seen today, from FOPT for BPF1, BPF2, BPF3 with
$v_{R}=10^{4},10^{5},10^{6}$ GeV respectively, in comparison to _power-law-
integrated-sensitivities_ [78] of different experiments. Figure 14: Strength
of GW peak for different values of $\alpha$ and $\beta/H_{*}$ for
$T_{n}=10^{4},10^{5},10^{6}$ GeV. Red dots signify unobservable peaks in BBO,
DECIGO, CE, ET, and LISA whereas blue dots implies observable signal.

We now summarise the results from all the above-listed contributions. Fig. 13
shows the GW spectrum from FOPT for the 3 benchmark points BPF1, BPF2, BPF3
with different $v_{R}$, as given in Table 2 along with the power-law-
integrated sensitivity curves of different gravitational wave experiments
mentioned in section 3.2. As discussed there, we see that the peak of BPF1 and
BPF2 are within the PLIS curve of BBO and DECIGO whereas BPF3 peak is outside
the PLIS curves of any of these experiments. Therefore we conclude that the
signals from BPF1 and BPF2 will be observable in BBO and DECIGO assuming a
threshold SNR of 1 in both, while BPF3 will not be observed. In Fig. 14, we
show GW peaks for typical values of $\alpha$ and $\beta/H_{*}$ obtained from
our model with different $T_{n}$, and color code them depending on whether
these peaks lie within the PLIS curve of any of these experiments. The red
dots correspond to observable peaks whereas the blue dots represent peaks
outside the range of BBO, DECIGO, CE, ET, and LISA. We see that for energy
$10^{6}$ GeV it is a lot difficult to get detectable peak, which is verified
by us while scanning the model parameter space to find benchmark points.

The FOPT bubble contributions need to be supplemented by the contribution from
the disintegration of the residual DW. For the benchmark points considered
here, it is possible to add this contribution from the CHL-SOPT calculation of
Fig. 2 since at least the mass scales considered there are the same as used
for the benchmark points here. In Fig. 15 we plot doubly peaked spectra from
typical strong FOPTs at different scales, for the benchmark points (BPDs) 1,
2, 3 given in Table 4. The two peaks correspond to the FOPT peak and the decay
of residual domain walls. In Fig. 16 we plot the same for a high scale BPD,
$\eta\sim 10^{10}$GeV. It is seen from Fig. 15 that a direct verification of
degenerate 2-field FOPT is marginally possible for low scale L-R models,
whereas a high scale case such as in Fig. 16 is beyond the reach of currently
planned experiments.

| BPD1 | BPD2 | BPD3 | BPD4
---|---|---|---|---
$T_{n}/$GeV | $10^{4}$ | $10^{5}$ | $10^{6}$ | $10^{10}$
$\rho_{1}$ | $0.1$ | $0.1$ | $0.1$ | $0.1$
$\mu_{3}^{2}\sim T_{n}^{2}\rho_{1}$ | $10^{7}$ | $10^{9}$ | $10^{11}$ | $10^{19}$
$\epsilon$ | $10^{-26}$ | $10^{-26}$ | $10^{-26}$ | $10^{-17}$
$\alpha$ | 0.1 | 0.1 | 0.1 | 0.1
$\beta/H_{*}$ | 1000 | 1000 | 1000 | 1000
Table 4: Benchmark points for typical strong FOPT. Here we assumed
$T_{n}\sim\eta$ hence from Eq. (2.11), $\mu_{3}^{2}\sim T_{n}^{2}\rho_{1}$.
The resulting doubly peaked spectra are shown in Fig.s 15 and 16. Figure 15:
Doubly peaked spectra from strong FOPT for benchmark points 1, 2, 3 given in
Table 4. The peak resulting from wall collision, sound waves, and MHD
turbulence lies generally within the sensitivity of BBO and DECIGO whereas the
peak from residual domain wall decay generically occurs at low frequencies
accessible to PTA experiments.  Figure 16: Similar to Fig. 15, for benchmark
point 4 of Table 4.

##  6 Dissociating Parity breaking from LR Symmetry Breaking

In the minimal model that we discussed so far, parity $\mathcal{P}$ and LR
symmetry are broken simultaneously where the right-handed triplet $\Delta_{R}$
takes a high VEV whereas $\Delta_{L}$ VEV is assumed to be negligible. This
assumption is based on the fact that Standard Model is predominantly left-
handed, and what breaks this parity is still unknown. To address this issue,
Chang et al [8, 7] described a way to dissociate parity from LR symmetry by
adding a scalar $\chi$ into the particle content, with negative parity, i.e.
$\chi\rightarrow-\chi$ under $D$-parity transformation. This $\chi$ adds extra
terms to our tree level potential,

$\displaystyle V_{\chi}$
$\displaystyle=-\mu_{\chi}^{2}\chi^{2}+\gamma_{1}\chi^{4}$ (6.1)
$\displaystyle V_{\chi\Delta}$ $\displaystyle=M\chi\text{
Tr}\left(\Delta_{L}^{\dagger}\Delta_{L}-\Delta_{R}^{\dagger}\Delta_{R}\right)+\gamma_{2}\chi^{2}\text{
Tr}\left(\Delta_{L}^{\dagger}\Delta_{L}+\Delta_{R}^{\dagger}\Delta_{R}\right)$
(6.2) $\displaystyle V_{\chi r}$ $\displaystyle=\gamma_{2}\chi^{2}\text{
Tr}\left(r^{\dagger}r\right)+\gamma_{2}^{1}\chi^{2}\left(\text{Det}r+\text{Det}r^{\dagger}\right)$
(6.3)

In this scenario $\chi$ takes a VEV, breaking parity at $M_{\text{parity}}$
scale before $SU(2)_{R}$ is broken at scale $M_{R}$. Depending on which minima
the field $\chi$ takes out of $\pm
M_{\text{parity}}=\mu_{\chi}/\sqrt{2\gamma_{1}}$ in a region, the Universe
becomes left or right-handed in that region, subject to minimal fine-tuning.
Since our Universe is left-handed, we want the right-handed regions/domains to
vanish eventually. This can be achieved again by making the domain walls
unstable by adding a small bias term into the $\chi$ potential. If
$M_{\text{parity}}\gg M_{R}$, the contribution from $SU(2)_{R}$ or $SU(2)_{L}$
breaking within left-like and right-like domains to the GW background will be
negligible compared to the contribution from the $\chi$ domain walls. Thus we
can assume that this type of dissociated high-energy parity breaking scheme
can also be responsible for detectable GW signals.

##  7 Conclusion

All our calculations have been carried out with a single representative field
as the argument of the effective potential. While we have paid attention to
the qualitatively distinct nature of the phase transitions, realistic results
need to be supplemented by additional factors. We have discussed the LRSM with
two $\Delta$ fields but our calculations focus on a single field for
simplicity of calculations. A supersymmetric version of the LRSM contains a
further doubling of these fields, in addition to a few others that might be
needed as in [109, 110, 43]. Thus the power radiated as GW gets a factor of at
least 2 if not more for any realistic model. However since the comparisons
with experiments are shown in logarithmic scales, we have fairly
representative results.

For the single field case, a final summary is provided by Fig. 17. It contains
the contour of peak frequency and peak amplitude in the parameter space of
$\eta$ and $\epsilon$. This shows the region of the parameter space where an
observable signal is produced.

Figure 17: Contour of peak frequency and amplitude in the parameter space of
$\eta$ and $\epsilon$. The red region corresponds to $\epsilon$ too large to
form domain walls with scaling evolution. Grey region corresponds to the case
where $\epsilon$ is small and domain walls dominate the energy density
eventually. Here we have taken $\rho_{1}=0.4$ and $g_{*}=134$.

In summary, we have emphasized the qualitatively different considerations
which arise for phase transitions when dealing with a theory with left-right
symmetry, and similar considerations will apply to the case of $SO(10)$ when
the pseudo-defect domain walls are taken into account [111]. In particular,
the SOPT case, usually not expected to leave behind any signature indeed
creates a GW signature, which should be accessible to experiments being
currently planned. If the parameter space signals an FOPT then in fact we get
two separate peaks, one from the usual bubble wall collision picture and an
additional peak at much lower frequencies arising from the residual DW that
arise due to the presence of two or more degenerate fields present in the
theory. This case is marginally accessible for low scale L-R models, Fig. 15.
A high scale case such as in Fig. 16 is beyond the reach of currently planned
experiments.

##  Appendix A Appendix

###  A.1 Rescaling the Field Equation

There’s only one energy scale in our theory, $\eta$. We can define all other
dimensional parameters in terms of $\eta$ to make our equation dimensionless,
eg: $t\rightarrow t\eta$, $x\rightarrow x\eta$, $r\rightarrow r\eta^{-1}$,
etc. We assume the evolution takes place entirely in the radiation dominated
era such that the scale factor and the Hubble parameter have time dependence

$\displaystyle a(t)$ $\displaystyle=a_{0}t^{1/2}$ (A.1) $\displaystyle H(t)$
$\displaystyle=\dot{a}/a=\frac{1}{2t}$ (A.2)

Thus, we make them dimensionless by $a\rightarrow a\eta^{1/2}$, $H\rightarrow
H\eta^{-1}$. The rescaled equation is obtained after dividing Eq. (4.1) by
$\eta^{3}$,

$\ddot{r}+3H(t)\dot{r}-\frac{\nabla^{2}}{a(t)^{2}}r+\rho_{1}r\left(r^{2}-1\right)+\epsilon\left(r^{2}-1\right)=0$
(A.3)

where all the terms represent their scaled versions now.

###  A.2 Calculation of Gravitational Waves

Let us look at the theory of GWs in an expanding universe. We use Green’s
method of [112, 113] in this section. We consider GWs to be sourced by tensor
perturbations in the spatially flat FLRW metric,

$ds^{2}=-dt^{2}+a^{2}(t)\left(\delta_{ij}+h_{ij}\right)dx^{i}dx^{j}$ (A.4)

where $a(t)$ is the scale factor, $h_{ij}$ is the tensor perturbation
satisfying the transverse-traceless (TT) condition $\partial_{i}h_{ij}=0$ and
which obeys the linearized Einstein equation,

$\ddot{h}_{ij}(t,\mathbf{x})+3H\dot{h}_{ij}(t,\mathbf{x})-\frac{\nabla^{2}}{a^{2}}h_{ij}(t,\mathbf{x})=\frac{16\pi
G}{a^{2}}T_{ij}^{\mathrm{TT}}(t,\mathbf{x})$ (A.5)

where dot denotes derivative with respect to cosmic time $t$ and $T^{TT}$ is
the TT part of the energy-momentum tensor. Our treatment is valid for any
arbitrarily expanding universe with,

$a(t)\propto\tau^{\alpha}\propto t^{\gamma}$ (A.6)

In the spatial Fourier space and in terms of conformal time $d\tau=dt/a(t)$,
this equation takes the form,

$h^{\prime\prime}_{ij}(t,\mathbf{k})+\frac{2\alpha}{\tau}h^{\prime}_{ij}(\tau,\mathbf{k})+k^{2}h_{ij}(\tau,\mathbf{k})=16\pi
GT^{TT}_{ij}(\tau,\mathbf{k})$ (A.7)

where prime notation denotes derivative with respect to conformal time. In
terms of the rescaled metric,

$\bar{h}_{ij}=ah_{ij}$ (A.8)

we obtain,

$\left[\frac{\partial^{2}}{\partial
x^{2}}+\left(1-\frac{4\nu^{2}-1}{4x^{2}}\right)\right]\bar{h}_{ij}(\tau,\mathbf{k})=\frac{16\pi
Ga(\tau)}{k^{2}}T^{TT}_{ij}(\tau,\mathbf{k})$ (A.9)

where $x=k\tau$, and $\nu$ is defined as,

$\nu=\alpha-\frac{1}{2}=\frac{3\gamma-1}{2(1-\gamma)}$ (A.10)

Our source $T^{TT}$ is active in the interval $\tau_{i}<\tau<\tau_{f}$. The
solution of Eq. (A.9) with initial condition
$\bar{h}_{ij}(\tau_{i})=\bar{h}^{\prime}_{ij}(\tau_{i})=0$ can be found using
Green’s method (for $\tau\leq\tau_{f}$),

$\bar{h}_{ij}(\tau,\mathbf{k})=\frac{8\pi^{2}G}{k^{2}}\int_{x_{i}}^{x}dy(yx)^{1/2}\left[N_{\nu}(x)J_{\nu}(y)-J_{\nu}(x)N_{\nu}(y)\right]a(y)T_{ij}^{\mathrm{TT}}(y,\mathbf{k})$
(A.11)

where $J_{\nu}(x)$ and $N_{\nu}(x)$ are the Bessel and Neumann functions. For
$\nu=1/2$, we get the solution for radiation dominated era.

When $\tau\geq\tau_{f}$, the source becomes negligible and the homogeneous
solution to Eq. (A.9) is a linear combination of $J_{\nu}(x)$ and
$N_{\nu}(x)$,

$\bar{h}_{ij}(\tau,\mathbf{k})=A_{ij}(\mathbf{k})(k\tau)^{1/2}J_{\nu}(k\tau)+B_{ij}(\mathbf{k})(k\tau)^{1/2}N_{\nu}(k\tau)\quad\left(\text{
for }\tau\geq\tau_{f}\right)$ (A.12)

Coefficients $A_{ij}(\mathbf{k}),B_{ij}(\mathbf{k})$ can be found by equating
Eq. (A.11) and (A.12) at $\tau=\tau_{f}$,

$\displaystyle A_{ij}(\mathbf{k})$
$\displaystyle=-\frac{8\pi^{2}G}{k^{2}}\int_{x_{i}}^{x_{f}}dx\sqrt{x}a(x)N_{\nu}(x)T_{ij}^{\mathrm{TT}}(x,\mathbf{k})$
$\displaystyle B_{ij}(\mathbf{k})$
$\displaystyle=\frac{8\pi^{2}G}{k^{2}}\int_{x_{i}}^{x_{f}}dx\sqrt{x}a(x)J_{\nu}(x)T_{ij}^{\mathrm{TT}}(x,\mathbf{k})$
(A.13)

Thus, if we can construct the $T^{TT}(x,\mathbf{k})$, we can calculate the
perturbation $h_{ij}(\tau,\mathbf{k})$. The energy density of GWs is given by
[114],

$\displaystyle\rho_{\mathrm{gw}}$ $\displaystyle=\frac{1}{32\pi
G}\left\langle\dot{h}_{ij}(t,\mathbf{x})\dot{h}_{ij}(t,\mathbf{x})\right\rangle$
$\displaystyle\simeq\frac{1}{32\pi
Ga^{4}(\tau)}\left\langle\bar{h}_{ij}^{\prime}(\tau,\mathbf{x})\bar{h}_{ij}^{\prime}(\tau,\mathbf{x})\right\rangle$
(A.14)

where we have ignored higher order terms of $aH$ for the second equality,
since we assume that the wavelength of GWs is smaller than the Hubble radius,
e.g. $k\tau\gg 1$.

A more useful quantity for observational cosmology is the fraction of energy
density of GWs to the total energy density of the universe at time $t$,
defined by,

$\Omega_{gw}(t)=\frac{1}{\rho_{c}(t)}\frac{d\rho_{gw}(t)}{d\ln{k}}$ (A.15)

where $\rho_{c}(t)$ is the critical density of the universe at time $t$. Note
that $k$ is the magnitude of the comoving wavevector. It is related to the
frequency observed at time $t$ as,

$f=\frac{k}{2\pi a(t)}$ (A.16)

where we set the scale factor at the initial time $t_{i}$ to be unity,
$a(t_{i})=1$.

####  A.2.1 Discretizing the Equation of Motion

We follow the notation and method of discretization from [115]. We also modify
the codes given in the reference to fit our requirements for the simulation.
We start with the equation of motion in Eq. (A.3), with boundary conditions,

$\displaystyle r(x,y,z,0)$ $\displaystyle=I(x,y,z)$ (A.17)
$\displaystyle\frac{\partial}{\partial t}r(x,y,z,0)$
$\displaystyle=V(x,y,z)=0$ (A.18) $\displaystyle r(0,y,z,t)$
$\displaystyle=r(L_{x},y,z,t)$ (A.19) $\displaystyle r(x,0,z,t)$
$\displaystyle=r(x,L_{y},z,t)$ (A.20) $\displaystyle r(x,y,0,t)$
$\displaystyle=r(x,y,L_{z},t)$ (A.21)

where $I(x,y,z)$ is the prescribed initial condition. $L_{x}=L_{y}=L_{z}=b$ is
the size of the comoving simulation box. The last three equations represent
periodic boundary conditions.

Our 3D mesh contains $N_{x}\times N_{y}\times N_{z}$ mesh points. We also
subdivide the time domain $[t_{ini},t_{fin}=T]$ into $N_{t}$ discrete
equidistant points $\left\\{t_{0}=T_{i}<t_{1}<t_{2}<...<T\right\\}$. Now the
field can be defined at each mesh point at a given time by the notation
$r(x_{i},y_{j},z_{k},t_{n})=r^{n}_{ijk}$, where superscript $n$ denotes the
time coordinate and subscripts $i,j,k$ denote $x,y,z$ coordinates
respectively.

To discretize the wave equation, we replace the derivatives with central
differences as follows,

$\displaystyle\frac{\partial}{\partial t}r^{n}_{ijk}$
$\displaystyle=[D_{2t}r]^{n}_{ijk}\approx\frac{r^{n+1}_{ijk}-r^{n-1}_{ijk}}{2\Delta
t}$ (A.22) $\displaystyle\frac{\partial^{2}}{\partial t^{2}}r^{n}_{ijk}$
$\displaystyle=[D_{t}D_{t}r]^{n}_{ijk}\approx\frac{r^{n+1}_{ijk}-2r^{n}_{ijk}+r^{n-1}_{ijk}}{\Delta
t^{2}}$ (A.23) $\displaystyle\frac{\partial^{2}}{\partial x^{2}}r^{n}_{ijk}$
$\displaystyle=[D_{x}D_{x}r]^{n}_{ijk}\approx\frac{r^{n}_{i+1jk}-2r^{n}_{ijk}+r^{n}_{i-1jk}}{\Delta
x^{2}}$ (A.24)

and similarly for $y$ and $z$ coordinates. $D_{t},D_{x}$ represents the
derivative operators in a compact notation. In this notation (omitting the
space-time indices), Eq. (A.3) becomes,

$D_{t}D_{t}r+3HD_{2t}r-\frac{1}{a^{2}}\left[D_{x}D_{x}+D_{y}D_{y}+D_{z}D_{z}\right]r+\rho_{1}r(r^{2}-1)+\epsilon\left(r^{2}-1\right)=0$
(A.25)

This, in discrete form, can be written as,

$\displaystyle\frac{r^{n+1}_{ijk}-2r^{n}_{ijk}+r^{n-1}_{ijk}}{\Delta
t^{2}}+3H^{n}\frac{r^{n+1}_{ijk}-2r^{n}_{ijk}+r^{n-1}_{ijk}}{\Delta
t^{2}}-\frac{1}{(a^{n})^{2}}\Bigg{[}\frac{r^{n}_{i+1jk}-2r^{n}_{ijk}+r^{n}_{i-1jk}}{\Delta
x^{2}}$ $\displaystyle+\frac{r^{n}_{ij+1k}-2r^{n}_{ijk}+r^{n}_{ij-1k}}{\Delta
y^{2}}+\frac{r^{n}_{ijk+1}-2r^{n}_{ijk}+r^{n}_{ijk-1}}{\Delta
z^{2}}\Bigg{]}+\rho_{1}r^{n}_{ijk}\left((r^{n}_{ijk})^{2}-1\right)$
$\displaystyle+\epsilon\left((r^{n}_{ijk})^{2}-1\right)=0$ (A.26)

Provided we know the state $r^{n}$ and $r^{n-1}$ at time slices $n$ and $n-1$,
we can solve this algebraic equation for $r^{n+1}$,

$\displaystyle r^{n+1}_{ijk}$
$\displaystyle=\left(1+\frac{3}{2}A^{n}\right)^{-1}\Bigg{[}\left\\{2-\rho_{1}\Delta
t^{2}\left((r^{n}_{ijk})^{2}-1\right)\right\\}r^{n}_{ijk}+\epsilon\Delta
t^{2}\left((r^{n}_{ijk})^{2}-1\right)$
$\displaystyle+\left(\frac{3}{2}A^{n}-1\right)r^{n-1}_{ijk}+B^{n}_{x}\left\\{r^{n}_{i+1jk}-2r^{n}_{ijk}+r^{n}_{i-1jk}\right\\}+B^{n}_{y}\left\\{...\right\\}+B^{n}_{z}\left\\{...\right\\}\Bigg{]}$
(A.27)

where $A^{n}=H^{n}\Delta t$, and $B^{n}_{x}=\frac{\Delta
t^{2}}{(a^{n})^{2}\Delta x^{2}}$ and similarly for $B^{n}_{y},B^{n}_{z}$.

Eq. (A.2.1) is not valid for $n=0$, which is the first time step in the
simulation. We note that the initial velocity of the field is given by,

$V^{0}_{ijk}=\frac{r^{1}_{ijk}-r^{-1}_{ijk}}{2\Delta t}\implies
r^{-1}_{ijk}=r^{1}_{ijk}-2\Delta tV^{0}_{ijk}$ (A.28)

Substituting this into Eq. (A.2.1), and solving for $r^{1}_{ijk}$ gives,

$\displaystyle r^{1}_{ijk}=\frac{1}{2}$
$\displaystyle\Big{[}\left\\{2-\rho_{1}\Delta
t^{2}\left({r^{0}}_{ijk}^{2}-1\right)\right\\}r^{0}_{ijk}+\epsilon\Delta
t^{2}\left({r^{n}}_{ijk}^{2}-1\right)+\left(\frac{3}{2}A^{0}-1\right)2\Delta
tV$
$\displaystyle+B^{0}_{x}\\{...\\}+B^{0}_{y}\\{...\\}+B^{0}_{z}\\{...\\}\Big{]}$
(A.29)

Eq. (A.2.1) along with (A.2.1) represent the recursive relations that we will
use in our simulation.

####  A.2.2 Periodic Boundary Condition

Periodic boundary condition is implemented on the 6 sides of the simulation
box. It is done to mimic the simulation over large distances and long times.
The wave exiting through one side of the simulation box enters the box
immediately through the opposite side, that way making the wave motion
periodic. The strategy to implement the periodic boundary condition is as
follows: apply open boundary condition on side A of the box, and equate the
opposite side B with A. That way, any displacement at side A will be instantly
reflected at side B. And we have 3 such AB pairs.

An open boundary condition means a wave can pass through the boundary without
any resistance. In 3D, the condition we implement is [115],

$D_{t}^{+}r+c_{x}D_{x}^{-}r+c_{y}D_{y}^{-}r+c_{z}D_{z}^{-}r=0\big{|}^{n}_{N_{x},N_{y},N_{z}}$
(A.30)

where $D$ represents the derivative operators and the superscript $\pm$
represents forward or backward finite difference. In our case,
$c_{x}=c_{y}=c_{z}=\tilde{c}=1/\sqrt{3}$ are wave velocities along $x,y,z$
directions. This equation is implemented at all times (denoted by superscript
$n$) at the 3 boundaries of the cube (denoted by subscript
$N_{x},N_{y},N_{z}$). Solving for $r^{n+1}_{N_{x},j,k}$ at $N_{x}$ we get,

$\displaystyle r^{n+1}_{N_{x},j,k}=$ $\displaystyle
r^{n}_{N_{x},j,k}-\frac{\bar{C}_{x}}{\sqrt{3}}\left(r^{n}_{N_{x},j,k}-r^{n}_{N_{x}-1,j,k}\right)-\frac{\bar{C}_{y}}{\sqrt{3}}\left(r^{n}_{N_{x},j,k}-r^{n}_{N_{x},j-1,k}\right)$
$\displaystyle-\frac{\bar{C}_{z}}{\sqrt{3}}\left(r^{n}_{N_{x},j,k}-r^{n}_{N_{x},j,k-1}\right)$
(A.31)

where $\bar{C}_{x,y,z}=\frac{\Delta t}{\Delta x,y,z}$ are the Courant numbers.
We get similar recursive relations for $r^{n+1}_{N_{y}/N_{z}}$ at
$N_{y}/N_{z}$. We use these equations to update the $N_{x},N_{y},N_{z}$
boundaries of the box while we equate the opposite sides, at
$-N_{x},-N_{y},-N_{z}$ with them.

####  A.2.3 Stability Condition

We can consider the coefficients $B_{x},B_{y},B_{z}$ from Eq. (A.2.1) to be
analogous to Courant numbers [115], and exert the stability condition for the
simulation to be,

$\max(B^{n}_{x})+\max(B^{n}_{y})+\max(B^{n}_{z})\leq 1$ (A.32)

As $B^{n}_{x}=\frac{\Delta t^{2}}{a^{n2}\Delta x^{2}}$, this gives,

$\Delta t\leq\min(a^{n})\left(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta
y^{2}}+\frac{1}{\Delta z^{2}}\right)^{-1/2}=\left(\frac{1}{\Delta
x^{2}}+\frac{1}{\Delta y^{2}}+\frac{1}{\Delta z^{2}}\right)^{-1/2}$ (A.33)

Simulation Procedure

  1. 1.

Define initial conditions as in Eq. (4.10) by discretizing it.

  2. 2.

Perform the first time step using Eq. (A.2.1).

  3. 3.

Update the field using Eq. (A.2.1) in iteration for $n>0$.

## Acknowledgments

We thank the authors of [104] for bringing recent developments to our
attention. ZAB is supported by an INSPIRE fellowship. The research of UAY is
supported by an Institute chair professorship.

## References

  * [1] B. P. Abbott and et. al., Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116 (2016), no. 6 61102\.
  * [2] A. Einstein, Näherungsweise Integration der Feldgleichungen der Gravitation, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin (1, 1916) 688–696.
  * [3] C. Caprini and D. G. Figueroa, Cosmological backgrounds of gravitational waves, Classical and Quantum Gravity 35 (8, 2018) 163001.
  * [4] B. P. Abbott and et. al., Exploring the sensitivity of next generation gravitational wave detectors, Classical and Quantum Gravity 34 (2, 2017) 044001.
  * [5] P. R. Saulson, Fundamentals of Interferometric Gravitational Wave Detectors.  World Scientific, 2nd. ed. ed., 2017.
  * [6] A. N. Lommen, Pulsar timing arrays: the promise of gravitational wave detection, Reports on Progress in Physics 78 (12, 2015) 124901\.
  * [7] D. Chang, R. N. Mohapatra, and M. K. Parida, Decoupling of Parity- and SU ( 2 ) R -Breaking Scales: A New Approach to Left-Right Symmetric Models, Physical Review Letters 52 (3, 1984) 1072–1075.
  * [8] D. Chang, R. N. Mohapatra, and M. K. Parida, New approach to left-right-symmetry breaking in unified gauge theories, Physical Review D 30 (9, 1984) 1052–1063.
  * [9] K. S. Babu and A. Patra, Higgs boson spectra in supersymmetric left-right models, Physical Review D 93 (3, 2016) 055030.
  * [10] L. Basso, B. Fuks, M. E. Krauss, and W. Porod, Doubly-charged Higgs and vacuum stability in left-right supersymmetry, Journal of High Energy Physics 2015 2015:7 2015 (7, 2015) 1–25.
  * [11] R. N. Mohapatra and J. C. Pati, "Natural" left-right symmetry, Physical Review D 11 (5, 1975) 2558–2561.
  * [12] G. Senjanovic and R. N. Mohapatra, Exact left-right symmetry and spontaneous violation of parity, Physical Review D 12 (9, 1975) 1502–1505.
  * [13] G. C. Branco and L. Lavoura, Natural CP breaking in left-right symmetric theories, Physics Letters B 165 (12, 1985) 327–332.
  * [14] J. Basecq, J. Liu, J. Milutinovic, and L. Wolfenstein, Spontaneous CP Violation in SU(2)-l X SU(2)-r X U(1)-($B^{-}$l) Models, Nucl. Phys. B 272 (1986) 145–157.
  * [15] J. F. Gunion, J. Grifols, A. Mendez, B. Kayser, and F. Olness, Higgs bosons in left-right-symmetric models, Physical Review D 40 (9, 1989) 1546–1561.
  * [16] G. Barenboim, M. Gorbahn, U. Nierste, and M. Raidal, Higgs Sector of the Minimal Left-Right Symmetric Model, Phys. Rev. D 65 (2002) 95003\.
  * [17] K. Kiers, M. Assis, and A. A. Petrov, Higgs sector of the left-right model with explicit C P violation, Physical Review D 71 (6, 2005) 115015.
  * [18] J. Chakrabortty, P. Konar, and T. Mondal, Copositive criteria and boundedness of the scalar potential, Physical Review D 89 (5, 2014) 095008.
  * [19] B. Brahmachari, M. K. Samal, and U. Sarkar, Potential Minimization in Left-Right Symmetric Models, .
  * [20] H. Fritzsch and P. Minkowski, Unified interactions of leptons and hadrons, Annals of Physics 93 (1975), no. 1 193–266.
  * [21] H. Georgi, The state of the art—gauge theories, in AIP Conference Proceedings, vol. 23, pp. 575–582, American Institute of Physics, 1975.
  * [22] C. Arbeláez, J. C. Romão, M. Hirsch, and M. Malinský, LHC-scale left-right symmetry and unification, Physical Review D 89 (2, 2014) 035002.
  * [23] F. F. Deppisch, T. E. Gonzalo, and L. Graf, Surveying the SO(10) model landscape: The left-right symmetric case, Physical Review D 96 (9, 2017) 055003.
  * [24] N. T. Shaban and W. J. Stirling, Minimal left-right symmetry and SO(10) grand unification using LEP coupling constant measurements, Physics Letters B 291 (9, 1992) 281–287.
  * [25] C.-Y. Chen, P. S. B. Dev, and R. N. Mohapatra, Probing heavy-light neutrino mixing in left-right seesaw models at the LHC, Physical Review D 88 (8, 2013) 033014.
  * [26] M. Lindner, F. S. Queiroz, W. Rodejohann, and C. E. Yaguna, Left-right symmetry and lepton number violation at the Large Hadron electron Collider, Journal of High Energy Physics 2016 (2016), no. 6 140.
  * [27] S. Patra, F. S. Queiroz, and W. Rodejohann, Stringent dilepton bounds on left-right models using LHC data, Physics Letters B 752 (1, 2016) 186–190.
  * [28] P. S. Dev, D. Kim, and R. N. Mohapatra, Disambiguating seesaw models using invariant mass variables at hadron colliders, Journal of High Energy Physics 2016 2016:1 2016 (1, 2016) 1–37.
  * [29] Y. Rodríguez and C. Quimbay, Spontaneous CP phases and flavour changing neutral currents in the left–right symmetric model, Nuclear Physics B 637 (8, 2002) 219–242.
  * [30] M. Nemevšek, F. Nesti, and G. Popara, Keung-Senjanović process at the LHC: From lepton number violation to displaced vertices to invisible decays, Physical Review D 97 (6, 2018) 115018.
  * [31] P. S. Dev, R. N. Mohapatra, and Y. Zhang, Probing the Higgs sector of the minimal Left-Right symmetric model at future hadron colliders, Journal of High Energy Physics 2016 2016:5 2016 (5, 2016) 1–55.
  * [32] A. Maiezza, M. Nemevšek, F. Nesti, and G. Senjanović, Left-right symmetry at LHC, Physical Review D 82 (9, 2010) 055022.
  * [33] R. Rajaraman, Solitons and instantons: An introduction to solitons and instantons in quantum field theory.-North-Holland Publishing Company: Amsterdam, New York, Oxford, .
  * [34] G. B. Gelmini, M. Gleiser, and E. W. Kolb, Cosmology of biased discrete symmetry breaking, Physical Review D 39 (3, 1989) 1558–1566.
  * [35] Y. B. Zel’dovich, I. Y. Kobzarev, and L. B. Okun, Cosmological consequences of spontaneous violation of discrete symmetry, Zh. Eksp. Teor. Fiz. 40 (1974) 3–11.
  * [36] B. Rai and G. Senjanović, Gravity and the domain-wall problem, Physical Review D 49 (3, 1994) 2729–2733.
  * [37] A. Vilenkin, Gravitational field of vacuum domain walls and strings, Physical Review D 23 (2, 1981) 852–857.
  * [38] P. Sikivie, Axions, Domain Walls, and the Early Universe, Physical Review Letters 48 (4, 1982) 1156–1159.
  * [39] A. K. Mohanty and F. W. Stecker, Matter-antimatter domains: A possible solution to the CP domain wall problem in the early universe, Physics Letters B 143 (8, 1984) 351–356.
  * [40] S. Mishra and U. A. Yajnik, Spontaneously broken parity and consistent cosmology with transitory domain walls, Phys. Rev. D 81 (2010) 045010, [arXiv:0911.1578].
  * [41] P. Banerjee and U. A. Yajnik, New ultraviolet operators in supersymmetric SO(10) GUT and consistent cosmology, Phys. Rev. D 101 (2020) 075041, [arXiv:1812.11475].
  * [42] J. M. Cline, U. A. Yajnik, S. N. Nayak, and M. Rabikumar, Transient domain walls and lepton asymmetry in the left-right symmetric model, Phys. Rev. D 66 (2002) 065001, [hep-ph/0204319].
  * [43] A. Sarkar, Abhishek, and U. A. Yajnik, PeV scale left-right symmetry and baryon asymmetry of the Universe, Nucl. Phys. B 800 (2008) 253–269, [arXiv:0710.5410].
  * [44] E. Witten, Cosmic separation of phases, Physical Review D 30 (7, 1984) 272–285.
  * [45] C. J. Hogan, Nucleation of cosmological phase transitions, Physics Letters B 133 (1983), no. 3 172–176.
  * [46] C. J. Hogan, Gravitational radiation from cosmological phase transitions, Monthly Notices of the Royal Astronomical Society 218 (2, 1986) 629–636.
  * [47] M. S. Turner and F. Wilczek, Relic gravitational waves and extended inflation, Physical Review Letters 65 (12, 1990) 3080–3083.
  * [48] M. Kamionkowski, A. Kosowsky, and M. S. Turner, Gravitational radiation from first-order phase transitions, Physical Review D 49 (3, 1994) 2837–2851.
  * [49] C. Grojean and G. Servant, Gravitational waves from phase transitions at the electroweak scale and beyond, Physical Review D 75 (2, 2007) 043507.
  * [50] J. Ellis, M. Lewicki, and J. M. No, On the maximal strength of a first-order electroweak phase transition and its gravitational wave signal, Journal of Cosmology and Astroparticle Physics 2019 (4, 2019) 003–003.
  * [51] V. Brdar, L. Graf, A. J. Helmboldt, and X.-J. Xu, Gravitational waves as a probe of left-right symmetry breaking, Journal of Cosmology and Astroparticle Physics 2019 (12, 2019) 027–027.
  * [52] D. Borah and A. Dasgupta, Probing left-right symmetry via gravitational waves from domain walls, Phys. Rev. D 106 (2022), no. 3 035016, [arXiv:2205.12220].
  * [53] M. Kawasaki and K. Saikawa, Study of gravitational radiation from cosmic domain walls, Journal of Cosmology and Astroparticle Physics 2011 (9, 2011) 008–008.
  * [54] T. Hiramatsu, M. Kawasaki, and K. Saikawa, On the estimation of gravitational wave spectrum from cosmic domain walls, .
  * [55] K. . Ichi Saikawa, A Review of Gravitational Waves from Cosmic Domain Walls, Universe 2017, Vol. 3, Page 40 3 (5, 2017) 40.
  * [56] R. N. Mohapatra and J. C. Pati, A Natural Left-Right Symmetry, Phys. Rev. D11 (1975) 2558.
  * [57] J. C. Pati and A. Salam, Lepton Number as the Fourth Color, Phys. Rev. D10 (1974) 275–289. [Erratum: Phys. Rev.D11,703(1975)].
  * [58] G. Senjanovic and R. N. Mohapatra, Exact Left-Right Symmetry and Spontaneous Violation of Parity, Phys.Rev. D12 (1975) 1502.
  * [59] G. Senjanovic, Spontaneous Breakdown of Parity in a Class of Gauge Theories, Nucl. Phys. B153 (1979) 334–364.
  * [60] R. N. Mohapatra and G. Senjanovic, Neutrino mass and spontaneous parity nonconservation, Phys. Rev. Lett. 44 (1980) 912.
  * [61] R. N. Mohapatra and G. Senjanovic, Neutrino Masses and Mixings in Gauge Models with Spontaneous Parity Violation, Phys. Rev. D23 (1981) 165.
  * [62] J. C. Pati and A. Salam, Unified Lepton-Hadron Symmetry and a Gauge Theory of the Basic Interactions, Phys. Rev. D8 (1973) 1240–1251.
  * [63] J. C. Pati and A. Salam, Are There Anomalous Lepton-Hadron Interactions?, Phys. Rev. Lett. 32 (1974) 1083.
  * [64] N. G. Deshpande, J. F. Gunion, B. Kayser, and F. Olness, Left-right-symmetric electroweak models with triplet Higgs field, Physical Review D 44 (8, 1991) 837–858.
  * [65] D. A. Kirzhnits and A. D. Linde, Macroscopic Consequences of the Weinberg Model, Phys. Lett. B 42 (1972) 471–474.
  * [66] S. Weinberg, Gauge and Global Symmetries at High Temperature, Phys. Rev. D 9 (1974) 3357–3378.
  * [67] L. Dolan and R. Jackiw, Symmetry Behavior at Finite Temperature, Phys. Rev. D 9 (1974) 3320–3341.
  * [68] T. W. B. Kibble, Topology of Cosmic Domains and Strings, J. Phys. A 9 (1976) 1387–1398.
  * [69] T. W. B. Kibble, Some Implications of a Cosmological Phase Transition, Phys. Rept. 67 (1980) 183.
  * [70] J. W. Essam, Percolation theory, Reports on Progress in Physics 43 (7, 1980) 833–912.
  * [71] D. Stauffer, Scaling theory of percolation clusters, Physics Reports 54 (7, 1979) 1–74.
  * [72] T. Hiramatsu, M. Kawasaki, and K. Saikawa, Gravitational waves from collapsing domain walls, Journal of Cosmology and Astroparticle Physics 2010 (5, 2010) 032–032.
  * [73] M. Li, Q.-S. Yan, Y. Zhang, and Z. Zhao, Prospects of gravitational waves in the minimal left-right symmetric model, Journal of High Energy Physics 2021 (2021), no. 3 267.
  * [74] C. L. Wainwright, CosmoTransitions: Computing cosmological phase transition temperatures and bubble profiles with multiple fields, Computer Physics Communications 183 (9, 2012) 2006–2013.
  * [75] K. Kadota, M. Kawasaki, and K. Saikawa, Gravitational waves from domain walls in the next-to-minimal supersymmetric standard model, Journal of Cosmology and Astroparticle Physics 2015 (10, 2015) 041–041.
  * [76] T. Hiramatsu, M. Kawasaki, and K. Saikawa, On the estimation of gravitational wave spectrum from cosmic domain walls, Journal of Cosmology and Astroparticle Physics 2014 (2, 2014) 031–031.
  * [77] A. Paul, U. Mukhopadhyay, and D. Majumdar, Gravitational wave signatures from domain wall and strong first-order phase transitions in a two complex scalar extension of the Standard Model, Journal of High Energy Physics 2021 2021:5 2021 (5, 2021) 1–28.
  * [78] K. Schmitz, New Sensitivity Curves for Gravitational-Wave Signals from Cosmological Phase Transitions, JHEP 01 (2021) 097, [arXiv:2002.04615].
  * [79] M. Maggiore, Gravitational wave experiments and early universe cosmology, Phys. Rept. 331 (2000) 283–367, [gr-qc/9909001].
  * [80] B. Allen, The Stochastic gravity wave background: Sources and detection, in Les Houches School of Physics: Astrophysical Sources of Gravitational Radiation, pp. 373–417, 4, 1996.  gr-qc/9604033.
  * [81] B. Allen and J. D. Romano, Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities, Phys. Rev. D 59 (1999) 102001, [gr-qc/9710117].
  * [82] W. H. Press, B. S. Ryden, and D. N. Spergel, Dynamical Evolution of Domain Walls in an Expanding Universe, Astrophys. J. 347 (1989) 590–604.
  * [83] D. Coulson, Z. Lalak, and B. A. Ovrut, Biased domain walls, Phys. Rev. D 53 (1996) 4237–4246.
  * [84] T. Garagounis and M. Hindmarsh, Scaling in numerical simulations of domain walls, Phys. Rev. D 68 (2003) 103506, [hep-ph/0212359].
  * [85] J. C. R. E. Oliveira, C. J. A. P. Martins, and P. P. Avelino, The Cosmological evolution of domain wall networks, Phys. Rev. D 71 (2005) 083509, [hep-ph/0410356].
  * [86] P. P. Avelino, J. C. R. E. Oliveira, and C. J. A. P. Martins, Understanding domain wall network evolution, Phys. Lett. B 610 (2005) 1–8, [hep-th/0503226].
  * [87] S. E. Larsson, S. Sarkar, and P. L. White, Evading the cosmological domain wall problem, Phys. Rev. D 55 (1997) 5129–5135, [hep-ph/9608319].
  * [88] M. Hindmarsh, Analytic scaling solutions for cosmic domain walls, Phys. Rev. Lett. 77 (1996) 4495–4498, [hep-ph/9605332].
  * [89] M. Hindmarsh, Level set method for the evolution of defect and brane networks, Phys. Rev. D 68 (2003) 043510, [hep-ph/0207267].
  * [90] P. P. Avelino, C. J. A. P. Martins, and J. C. R. E. Oliveira, One-scale model for domain wall network evolution, Phys. Rev. D 72 (2005) 083506, [hep-ph/0507272].
  * [91] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Gravitational Waves from the Sound of a First Order Phase Transition, Physical Review Letters 112 (1, 2014) 041301.
  * [92] S. Coleman, Fate of the false vacuum: Semiclassical theory, Physical Review D 15 (5, 1977) 2929–2936.
  * [93] C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, A. Petiteau, P. Schwaller, G. Servant, and D. J. Weir, Science with the space-based interferometer eLISA. II: gravitational waves from cosmological phase transitions, Journal of Cosmology and Astroparticle Physics 2016 (4, 2016) 001–001.
  * [94] M. Lewicki and V. Vaskonen, Gravitational waves from bubble collisions and fluid motion in strongly supercooled phase transitions, Eur. Phys. J. C 83 (2023), no. 2 109, [arXiv:2208.11697].
  * [95] J. Ellis, M. Lewicki, and V. Vaskonen, Updated predictions for gravitational waves produced in a strongly supercooled phase transition, JCAP 11 (2020) 020, [arXiv:2007.15586].
  * [96] M. Kierkla, A. Karam, and B. Swiezewska, Conformal model for gravitational waves and dark matter: a status update, JHEP 03 (2023) 007, [arXiv:2210.07075].
  * [97] A. D. Linde, Decay of the false vacuum at finite temperature, Nuclear Physics B 216 (5, 1983) 421–445.
  * [98] J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, Energy budget of cosmological first-order phase transitions, Journal of Cosmology and Astroparticle Physics 2010 (6, 2010) 028–028.
  * [99] A. H. Guth and E. J. Weinberg, Cosmological consequences of a first-order phase transition in the s${\mathrm{u}}_{5}$ grand unified model, Phys. Rev. D 23 (Feb, 1981) 876–885.
  * [100] V. K. S. Shante and S. Kirkpatrick, An introduction to percolation theory, Advances in Physics 20 (1971) 325–357.
  * [101] A. Kosowsky, M. S. Turner, and R. Watkins, Gravitational radiation from colliding vacuum bubbles, Physical Review D 45 (6, 1992) 4514–4535.
  * [102] C. Caprini and R. Durrer, Gravitational waves from stochastic relativistic sources: Primordial turbulence and magnetic fields, Physical Review D 74 (9, 2006) 063521.
  * [103] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Numerical simulations of acoustically generated gravitational waves at a first order phase transition, Physical Review D 92 (12, 2015) 123009.
  * [104] H.-K. Guo, K. Sinha, D. Vagie, and G. White, Phase transitions in an expanding universe: stochastic gravitational waves in standard and non-standard histories, Journal of Cosmology and Astroparticle Physics 2021 (jan, 2021) 001.
  * [105] M. Lewicki, M. Merchand, and M. Zych, Electroweak bubble wall expansion: gravitational waves and baryogenesis in Standard Model-like thermal plasma, JHEP 02 (2022) 017, [arXiv:2111.02393].
  * [106] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Shape of the acoustic gravitational wave power spectrum from a first order phase transition, Phys. Rev. D 96 (2017), no. 10 103520, [arXiv:1704.05871]. [Erratum: Phys.Rev.D 101, 089902 (2020)].
  * [107] C. Caprini, R. Durrer, and G. Servant, The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition, Journal of Cosmology and Astroparticle Physics 2009 (12, 2009) 024–024.
  * [108] P. Binétruy, A. Bohé, C. Caprini, and J.-F. Dufaux, Cosmological backgrounds of gravitational waves and eLISA/NGO: phase transitions, cosmic strings and other sources, Journal of Cosmology and Astroparticle Physics 2012 (6, 2012) 027–027.
  * [109] C. S. Aulakh, A. Melfo, and G. Senjanovic, Minimal supersymmetric left-right model, Phys. Rev. D57 (1998) 4174–4178. hep-ph/9707256.
  * [110] C. S. Aulakh, K. Benakli, and G. Senjanovic, Reconciling supersymmetry and left-right symmetry, Phys. Rev. Lett. 79 (1997) 2188–2191. hep-ph/9703434.
  * [111] I. Garg and U. A. Yajnik, Topological pseudodefects of a supersymmetric $SO(10)$ model and cosmology, Phys. Rev. D 98 (2018), no. 6 063523, [arXiv:1802.03915].
  * [112] D. Stauffer, Scaling theory of percolation clusters, Phys. Rept. 54 (1979) 1–74.
  * [113] J.-F. Dufaux, G. Felder, L. Kofman, and O. Navros, Gravity Waves from Tachyonic Preheating after Hybrid Inflation, JCAP 03 (2009) 001, [arXiv:0812.2917].
  * [114] M. Maggiore, Gravitational waves: Volume 1: Theory and experiments, vol. 1.  Oxford university press, 2008.
  * [115] H. P. Langtangen, Finite difference methods for wave motion, Department of Informatics, University of Oslo, preliminary version edition (2016).