††institutetext: Institute of Theoretical Physics, Faculty of Physics,
University of Warsaw, ul. Pasteura 5,
PL-02-093 Warsaw, Poland

# First-order phase transitions in Twin Higgs models

Marcin Badziak  and Ignacy Nałęcz<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We revisit phase transitions in Twin Higgs (TH) models. We show that strong
first-order phase transitions (FOPTs) can occur provided that appropriate
source of $\mathbb{Z}_{2}$ symmetry breaking between the twin and Standard
Model (SM) sectors is present. We found FOPTs in two classes of models. First:
with hard $\mathbb{Z}_{2}$ breaking in the scalar potential allowing for FOPT.
Second: with $\mathbb{Z}_{2}$ broken by enhanced Yukawa couplings of twin
leptons. We also considered supersymmetric UV completion of the second
scenario with light sleptons. The signal of gravitational waves produced
during these phase transitions is typically small but can be close to the
reach of AEDGE and Einstein Telescope in the case of the FOPT induced by light
twin sleptons. Our results open a way to generate SM baryon asymmetry in TH
models.

##  1 Introduction

Twin Higgs (TH) models Chacko:2005pe ; Chacko:2005vw ; Chacko:2005un
elegantly solve the little hierarchy problem of the SM by introducing a twin
(or mirror) copy of the Standard Model (SM). In this framework, the SM-like
Higgs boson is a pseudo-Nambu-Goldstone boson (PNGB) of a global SU(4)
symmetry and the $\mathbb{Z}_{2}$ symmetry relating the SM and its twin
counterparts allows to avoid breaking of the SU(4) symmetry by quantum
corrections. Importantly, the top partners responsible for the naturalness of
the electroweak (EW) scale are charged under twin SU(3) gauge group rather
than the SM color gauge group, easily avoiding bounds from the LHC. In the UV
completions of the TH model utilizing supersymmetry (SUSY) Falkowski:2006qq ;
Chang:2006ra ; Craig:2013fga ; Katz:2016wtw ; badziak2017S ; badziak2017M ;
badziak2017A , higher dimensions Craig:2014aea ; Craig:2014roa ;
Geller:2014kta  or composite Higgs Batra:2008jy ; Barbieri:2015lqa ;
Low:2015nqa ; Csaki:2015gfd ; Contino:2017moj  to solve the big hierarchy
problem of the SM, the mass scale of colored top partners can be pushed up to
a few-TeV range waiting for the discovery in a future collider.

Even though colored top partners in TH models are expected to be beyond the
reach of the LHC, this class of models can be tested at the LHC via the Higgs
coupling measurements and searches for additional Higgs bosons Barbieri:2005ri
; Craig:2015pha ; Buttazzo:2015bka ; Katz:2016wtw ; Ahmed:2017psb . Even more
signatures of TH models are expected in cosmological context due to the
richness of the twin sector. Several good dark matter (DM) candidates residing
in the twin sector have been identified that can be searched for in direct-
detection experiments Barbieri:2005ri ; Craig:2015pha ; GarciaGarcia:2015fol ;
Farina:2015uea ; Craig:2015xla ; Barbieri:2016zxn ; Barbieri:2017opf ;
Badziak:2019zys ; Terning:2019hgj ; Chacko:2021vin ; Badziak:2022eag . Mirror
sector of TH models may also modify the cosmic microwave background (CMB) and
the large-scale structure Prilepina:2016rlq ; Chacko:2018vss ; Bansal:2021dfh
. Moreover, the twin sector also generically contributes to dark radiation. As
a matter of fact, this contribution is often too big to be consistent with the
data and additional mechanism is required to solve this dark radiation
problem. Several solutions to this problem have been proposed, including late
entropy production that dilutes dark radiation Chacko:2016hvu ; Craig:2016lyx
, Fraternal Twin Higgs Craig:2015pha  where the first two generations of
fermions are removed or $\mathbb{Z}_{2}$ breaking in the Yukawa couplings of
light fermions Barbieri:2016zxn ; Barbieri:2017opf ; Harigaya:2019shz .

In the present article, we investigate phase transitions in TH models and ask
a question whether any of such transitions can be first-order. It is well
known that first-order phase transition (FOPT) may result in production of
gravitational waves (GW) Kamionkowski:1993fg ; Ellis:2018mja ; Caprini:2015zlo
; Caprini:2019egz  and may also allow for baryogenesis Kuzmin:1985mm ;
Shaposhnikov:1986jp , for a review see e.g. Cline:2006ts ; Morrissey:2012db .
Therefore, a presence of strong FOPT in the cosmological evolution of TH
models would provide another handle to test this class of SM extensions.

The EW phase transition in the SM is known to be smooth Arnold:1992rz ;
Kajantie:1995kf ; Kajantie:1996mn ; Kajantie:1996qd  but extending the Higgs
sector by adding a singlet may result in the first-order EW phase transition,
see e.g. refs.Barger:2007im ; Brauner:2016fla ; Curtin:2014jma . The Higgs
sector of TH models consists of the twin Higgs boson which is a singlet under
the SM gauge group. However, the structure of the scalar potential is highly
restricted by the approximate SU(4) and $\mathbb{Z}_{2}$ symmetries so it is
not obvious whether such extension allows for the first-order phase
transition.

Phase transitions in TH models have been studied in ref. Fujikura:2018duw
reaching the conclusion that they are not first-order. However, the analysis
of ref. Fujikura:2018duw  assumed that twin EW phase transition, in which twin
Higgs field obtained vacuum expectation value (VEV), occurred first and it
finished before the start of the EW transition. Such assumption allowed for
investigating twin EW and standard EW phase transitions separately and only
one field was dynamical at a time.

In the present paper we show that one-field approximation for phase
transitions is not always valid and studying two-field evolution is necessary
to get reliable predictions. We also find that strong FOPT can be present if
one includes effects of $\mathbb{Z}_{2}$ symmetry breaking. We demonstrate
this in two cases. One: when hard $\mathbb{Z}_{2}$ breaking is present in the
tree-level potential. Two: when $\mathbb{Z}_{2}$ symmetry is broken in Yukawa
couplings of light fermions.

In the first scenario, we revisited the parameter space of the TH model with
an unspecified UV completion (we refer to this scenario as UV agnostic TH) and
found a new region where the naturalness can be combined with viable SM Higgs
phenomenology, which was not considered in many previous analyzes Katz:2016wtw
; Fujikura:2018duw . The FOPT can be present only in this new region which is
characterized by hard $\mathbb{Z}_{2}$ breaking and negative
$\mathbb{Z}_{2}$-preserving but $SU(4)$-breaking quartic coupling.

In the second scenario, we found FOPTs induced by heavy twin leptons for some
range of twin lepton Yukawa couplings. In a supersymmetric UV completion of
this scenario, we found that light sleptons can induce even stronger FOPT. In
some part of parameter space the EW symmetry is not restored for temperatures
up to above 1 TeV. The presence of FOPT and EW symmetry non-restoration up to
high temperatures, may have important implications for baryogenesis.

We also present numerically simulated spectra of stochastic GW background
which can be emitted in the FOPT induced either by hard
$\mathbb{Z}_{2}$-breaking in the scalar potential or by asymmetric Yukawa
couplings. In the UV agnostic TH, we found the GW signal much too weak to be
detected in the foreseeable future. On the other hand, in supersymmetric
extension of TH with light sleptons and enhanced Yukawa couplings of twin
leptons the GW signal is much stronger but may be within reach of the future
GW detectors only in a small corner of the parameter space.

The rest of the article is organized as follows. In section 2 we introduce the
TH model and its free parameter space, and recall the one-loop effective
potential together with leading-order thermal corrections to the scalar
potential. Section 3 contains discussion of two possible scenarios of the EW
and twin EW symmetry breaking and their dependence on explicit
$\mathbb{Z}_{2}$ symmetry breaking between the SM and twin sector. In section
4 we present the numerical evidence for the strong FOPT taking place in SM and
twin sectors in UV agnostic TH model with hard $\mathbb{Z}_{2}$ breaking and
identify regions in the parameter space where such transition occurs. Section
5 is devoted to the analysis of phase transitions in the UV agnostic TH model
with $\mathbb{Z}_{2}$ breaking in lepton Yukawa couplings, whereas in section
6 we consider supersymmetric TH model with light sleptons and large twin
lepton Yukawa couplings and analyze their impact on scalar phase transitions.
Section 7 covers the computation of the GW signal emitted in the FOPTs that
were found in the preceding sections. Finally, section 8 contains the
discussion and summary of our results, as well as some comments on possible
future work.

##  2 Twin Higgs scalar potential

###  2.1 The tree-level potential

The model contains two scalar doublets: $H_{A}$ which transforms under the SM
$SU(2)_{L}\times U(1)_{Y}$ gauge group and its twin $H_{B}$ which transforms
under the twin gauge group $SU(2)_{L^{\prime}}\times U(1)_{Y^{\prime}}$. The
most general tree-level potential of the TH model reads

$V_{\text{tree}}(H_{A},H_{B})=\lambda(|H_{A}|^{2}+|H_{B}|^{2}-f_{0}^{2}/2)^{2}+\kappa(|H_{A}|^{4}+|H_{B}|^{4})+f_{0}^{2}\sigma|H_{A}|^{2}+\rho|H_{A}|^{4}\,,$
(1)

with

$H_{I}=\begin{pmatrix}\chi^{+}_{I}\\\
\frac{1}{\sqrt{2}}(h_{I}+i\,\chi^{0}_{I})\end{pmatrix}\,.$

Further on we exploit the unitary gauge where $\chi^{+}_{I}$ and
$\chi^{0}_{I}$ are eaten by longitudinal degrees of freedom of the SM and twin
$W^{\pm}$ and $Z^{0}$ vectors. Thus, one can substitute
$H_{I}=\left(\begin{matrix}0\\\ \frac{h_{I}}{\sqrt{2}}\end{matrix}\right)$
into (1), obtaining

$V_{\text{tree}}(h_{A},h_{B})=\frac{\lambda}{4}(h_{A}^{2}+h_{B}^{2}-f_{0}^{2})^{2}+\frac{\kappa}{4}(h_{A}^{4}+h_{B}^{4})+\frac{f_{0}^{2}\sigma}{2}h_{A}^{2}+\frac{\rho}{4}h_{A}^{4}\,.$
(2)

Here, $\lambda$, $\kappa$, $\sigma$, $\rho$ and $f_{0}$ are free real
parameters. $f_{0}$ is the scale of the spontaneous $SU(4)$ symmetry breaking
and is related to the VEVs of SM and twin Higgs bosons: $f_{0}^{2}\approx
v_{A}^{2}+v_{B}^{2}$, where $v_{I}\equiv\langle h_{I}\rangle$. While $f_{0}$
and $\lambda$ must be positive, there are generically no restriction on the
sign of the couplings that explicitly break the $SU(4)$ symmetry 111In the
past literature an arbitrary sign of $\sigma$ was considered while $\kappa$
was assumed to be positive Barbieri:2005ri ; Beauchesne:2015lva ; Katz:2016wtw
. This is well justified as long as one considers hard $\mathbb{Z}_{2}$
breaking coupling $\rho\leq 0$. However, in the general case there is no
reason to exclude negative values of $\kappa$ (see for example ref.
Ahmed:2017psb ).. In order to guarantee the approximate $SU(4)$ and
$\mathbb{Z}_{2}$ symmetries which protect the SM Higgs mass parameter from the
quadratically divergent loop corrections, the first term in the potential (2)
has to be dominant. Hence, the coupling $\lambda$ should dominate over the
other dimensionless scalar couplings: $|\kappa|$, $|\sigma|$, and $|\rho|$.
Nevertheless, the symmetry breaking interactions in the tree-level potential
play crucial roles: the $SU(4)$ breaking but $\mathbb{Z}_{2}$ preserving
$\kappa$ provides the non-zero mass of the SM-like Higgs while the softly
$\mathbb{Z}_{2}$-breaking $\sigma$ is necessary to obtain a mild hierarchy
$v_{B}\gtrsim 3v_{A}$ between the VEVs of the SM and twin Higgs bosons, which
is required by the LHC Higgs coupling measurements Ahmed:2017psb . The hard
$\mathbb{Z}_{2}$ breaking interaction controlled by $\rho$ is not absolutely
necessary and was often neglected in the past but it can be present in some UV
completions of TH models Katz:2016wtw .

The full expressions for the tree-level VEVs and masses of SM and twin Higgs
bosons were given in appendix A. Here, for better clarity, we present the
formulas obtained in the limit $\lambda\gg\kappa,\;\sigma,\;\rho$ and
$v_{A}/f_{0}\ll 1$ which is often refereed to as the pseudo Nambu-Goldstone
boson approximation. They read

$v_{A}^{2}\approx f_{0}^{2}\frac{\kappa-\sigma}{4\kappa+2\rho},\qquad
v_{B}^{2}\approx f_{0}^{2}-v_{A}^{2}\,,$ (3)

and

$m_{h}^{2}\approx 2v_{A}^{2}(2\kappa+\rho),\qquad m_{h^{\prime}}^{2}\approx
2f_{0}^{2}\lambda,$ (4)

where $m_{h}$ denotes SM-like Higgs boson mass, and $m_{h^{\prime}}$
corresponds to its twin counterpart.

By setting $v_{A}=246\;\text{GeV}$ and $m_{h}=125\;\text{GeV}$ one fixes tree-
level values of two out of five free parameters 222The elimination of free
parameters with the one-loop conditions for SM Higgs mass and VEV at $T=0$ is
a tedious numerical task. In spite of this, one-loop improved conditions were
used for all presented numerical results.. It is convenient to eliminate
$\sigma$ and $\rho$ and work with $\lambda$, $\kappa$ and $f_{0}$.

The $SU(4)$-breaking scale $f_{0}$ lays in a very specific range. The
consistency with the LHC data on the production and decays of the 125 GeV
Higgs requires $f_{0}\gtrsim 3v_{A}$ Ahmed:2017psb . On the other hand, in TH
models tuning grows with $f_{0}^{2}/v_{A}^{2}$ so to preserve naturalness of
the model $f_{0}$ should not be too large. For example, avoiding tuning worse
than 10 % requires $f_{0}\lesssim 5v_{A}$.

###  2.2 Loop and thermal corrections

In order to trace phase transition path, one should compute full effective
potential, consisting of zero temperature loop and thermal corrections to the
tree-level part (2)

$V_{\text{tot}}(h_{A},\;h_{B})=V_{\text{tree}}(h_{A},\;h_{B})+V_{\text{therm}}(h_{A},\;h_{B})+V_{\text{CW}}(h_{A},\;h_{B}).$
(5)

Throughout this work we use the Colman-Weinberg one-loop approximation for
zero temperature corrections, computed in the $\overline{\text{MS}}$
renormalization scheme Coleman:1973jx

$\displaystyle
V_{\text{CW}}(h_{A},\;h_{B})=\sum_{i\in\text{bosons}}\frac{n_{i}}{64\pi^{2}}m^{4}_{i}(h_{A},\;h_{B})\left(\log(\frac{m^{2}_{i}(h_{A},\;h_{B})}{\mu^{2}})-c_{i}\right)$
(6)
$\displaystyle-\sum_{i\in\text{fermions}}\frac{n_{i}}{64\pi^{2}}m^{4}_{i}(h_{A},\;h_{B})\left(\log\frac{m^{2}_{i}(h_{A},\;h_{B})}{\mu^{2}}-\frac{3}{2}\right),$

where $n$ is the number of the field degrees of freedom and the coefficients
$c_{i}$ are equal to $5/6$ in case of vectors and $3/2$ in case of scalars. We
set the renormalization scale to $\mu=f_{0}$ and computed the one loop
corrections to tree-level values of $\kappa$ and $\sigma$ in order to preserve
correct mass and the VEV of the SM Higgs.

The one-loop effective potential (6) and all the observables derived from it
depend on the choice of renormalization scale $\mu$. This is unfortunately
unavoidable if one works with the leading order approximation of the effective
potential. To estimate the renormalization scale dependence of our results we
computed the temperature at which global minimum becomes metastable, called
critical temperature $T_{c}$, using different choices of $\mu$. The critical
temperature is the physically meaningful quantity, since at this temperature
one expects the onset of the scalar phase transition in the universe. In the
FOPT that we discuss in section 3 for $\mu\in[200\text{ GeV},\,2000\text{
GeV}]$ the critical temperature computed using loop and thermally corrected
potential (5) varies by no more than $15\%$. This strongly indicates that any
reasonable choice of renormalization scheme leads to roughly the same physical
insight.

Thermal corrections to the scalar potential are given by Quiros:1999jp

$V_{\text{therm}}(h_{A},\;h_{B})=\frac{T^{4}}{2\pi^{2}}\sum_{i\in\text{bosons}}n_{i}\text{J}_{B}(\frac{m^{2}_{i}}{T^{2}})-\frac{T^{4}}{2\pi^{2}}\sum_{i\in\text{fermions}}n_{i}\text{J}_{F}(\frac{m^{2}_{i}}{T^{2}})$
(7)

where

$\displaystyle
J_{B}(x)=\int_{0}^{\infty}dk\,k^{2}\log[1-\exp(-\sqrt{k^{2}+x})],$ (8)
$\displaystyle
J_{F}(x)=\int_{0}^{\infty}dk\,k^{2}\log[1+\exp(-\sqrt{k^{2}+x})].$ (9)

At high temperatures the effective potential receives large contributions from
multi-loop bosonic diagrams which diverge in the infrared. This problem is
usually solved by the resummation of the leading-order contributions of these
diagrams Arnold:1992rz

$V_{\text{ring}}=-\sum_{i\in\text{bosons}}\frac{n_{i}T}{12\pi}[(\overline{m}^{2}_{i})^{\frac{3}{2}}-(m^{2}_{i})^{\frac{3}{2}}],$
(10)

or, more exactly, by replacing tree level masses in $V_{\text{therm}}$ and
$V_{\text{CW}}$ by their thermal analogues $\overline{m}^{2}_{i}$
Parwani:1991gq . In this work, we adopt the latter procedure. However, even
after the resummation we expect our one-loop approximation of the effective
potential to be robust if

$\frac{h_{I}(T)}{T}>g.$ (11)

Here, $h_{I}(T)$ denotes the Higgs VEV in the sector with broken
$\text{SU}(2)$ symmetry at temperature $T$ and $g$ is the weak coupling. For
more extensive discussion of the resummation and criterion (11) check appendix
B.2.

##  3 Phase dynamics in Twin Higgs models

The classification of possible scenarios for the evolution of fields is an
essential step towards finding the regions in parameter space, where the
scalar transition leads to events interesting from phenomenological point of
view. In general, during thermal evolution of the Universe the expectation
values of dynamical scalars follow some specific path in $h_{A}-h_{B}$ plane.
The scalar phase transition order and dynamics critically depend on its shape.

We assume that in the considered model at sufficiently high temperatures the
$SU(4)$ symmetry is restored, i.e. the global potential minimum is at
$h_{A}=h_{B}=0$ 333As far as the one-loop approximation of the effective
potential is considered, one could show that at sufficiently high temperatures
thermal corrections to the scalar masses would dominate field-dependent
contribution and force the symmetry restoration Kilic:2015joa .. Our aim is
now to determine which of the fields $h_{A}$ or $h_{B}$ acquires its
expectation value first. It is a well-known fact that in the non-abelian gauge
field theories the perturbative expansion of the effective potential breaks
down at high temperatures at which the symmetry is restored Arnold:1994bp ;
Linde:1980ts . The major problem is that we do not control the radiative
corrections from the transverse modes of vector fields. Thus, the perturbative
approximation of the effective potential (5) is, in general, inappropriate
tool for analyzing phase transitions at the origin.

Nevertheless, a direction at which the $\text{SU}(4)$ symmetry is broken
depends solely on the $\mathbb{Z}_{2}$-breaking part of the effective
potential while the large corrections from heavy vectors are
$\mathbb{Z}_{2}$-symmetric444This is true only if the SM and TS gauge
couplings are identical. Apparently, the TH naturalness requires them to be
very close, so with no lost of generality, one can set them formally equal..
For this reason, in the symmetric phase and its vicinity, one can simply
ignore large, unknown contributions from gauge sectors and use the one-loop
resummed effective potential (5) which encodes all the necessary information
about the $\mathbb{Z}_{2}$-breaking in the Lagrangian.

At the temperatures at which symmetry is restored, the effective potential is
well approximated by the polynomial in fields (see appendix B.1). The cubic
terms will be tiny, because at high temperatures mass of each particle is
dominated by its thermal corrections. Without significant cubic terms, the
field cannot tunnel from symmetric phase, and the only way to escape is the
continues phase transition.

When the continues transition takes place, the global minimum moves away from
the origin. This can happen only if the second order derivative computed at
the origin changes sign from positive to negative

$\displaystyle\frac{\partial^{2}V_{\text{tot}}}{\partial
h_{A}^{2}}\bigg{|}_{h_{A}=0,\,h_{B}=0}=\>\zeta_{A}T^{2}-\lambda
f_{0}^{2}(1-\frac{\sigma}{\lambda})$ (12)
$\displaystyle\frac{\partial^{2}V_{\text{tot}}}{\partial
h_{B}^{2}}\bigg{|}_{h_{A}=0,\,h_{B}=0}=\>\zeta_{B}T^{2}-\lambda f_{0}^{2}$
(13)

where $\zeta_{A}$ and $\zeta_{B}$ are thermal mass coefficients which depend
on the SM and TS couplings, respectively. The mixed derivative obviously
vanishes.

As the Universe cools down, both derivatives decrease and eventually one of
them reaches zero. At this point, continues phase transition occurs, and the
minimum leaves the origin and starts moving in $h_{A}$-$h_{B}$ plane. Let
$T_{A}$ and $T_{B}$ denote temperatures at which mass of the SM-like and twin
Higgs boson, respectively, would vanish. Then, one may compute the ratio

$\frac{T_{A}}{T_{B}}=\sqrt{\frac{\zeta_{B}}{\zeta_{A}}}\sqrt{1-\frac{\sigma}{\lambda}}.$
(14)

By evaluating the right-hand side of the above equation one finds which of the
fields $h_{A}$ or $h_{B}$ will first acquire non-zero VEV. The above reasoning
was presented in ref. Fujikura:2018duw . However, ref. Fujikura:2018duw
considered only the case in which $T_{B}\gg T_{A}$. eq. (14) shows that
$T_{B}>T_{A}$ if $\sigma>0$ and $\zeta_{A}\sim\zeta_{B}$ i.e. neglecting
effects of $\mathbb{Z}_{2}$ breaking in thermal masses. In our analysis we
relax the assumption $T_{B}\gg T_{A}$ and show that in the presence of
$\mathbb{Z}_{2}$ breaking an analysis of the two-field dynamics is necessary
to understand phase transitions in TH models.

It may seem that the order in which scalar fields acquire expectation values
is irrelevant since very soon the VEV of the remaining scalar will also start
growing. This, however, seems not to be the case. Once minimum leaves the
origin, the equations (12) and (13) should be corrected (note, that the mixed
derivatives are still zero)555Here, one can safely neglect zero temperature
loop corrections, since their contribution is sub-leading close to the
temperature of $SU(4)$-breaking transition.

$\displaystyle\frac{\partial^{2}V_{\text{tot}}}{\partial
h_{A}^{2}}\bigg{|}_{h_{A}=0,\,h_{B}}=\>\tilde{\zeta}_{A}T^{2}-\lambda
f_{0}^{2}(1-\frac{\sigma}{\lambda}-\frac{h_{B}^{2}}{f_{0}^{2}}),$ (15)
$\displaystyle\frac{\partial^{2}V_{\text{tot}}}{\partial
h_{B}^{2}}\bigg{|}_{h_{A},\,h_{B}=0}=\>\tilde{\zeta}_{B}T^{2}-\lambda
f_{0}^{2}(1-\frac{h_{A}^{2}}{f_{0}^{2}}).$ (16)

From the above equations it follows that non-zero VEV of one field stabilizes
the other one by giving positive contribution to its effective mass. Hence,
close to the origin the yet unbroken EW or twin EW symmetry can be broken only
by the first-order phase transition. On the other hand, when the difference
$|h_{A}-h_{B}|$ becomes large, thermal loop contributions from transverse
modes are not $\mathbb{Z}_{2}$-symmetric anymore. Thus, for big scalar
expectation values the above equation does not hold, and continues phase
transition is possible. In either case, critical temperature, transition order
and distance between degenerate minima strongly depend on the sector in which
the $\text{SU}(2)$ symmetry is broken first. This in turn is controlled by the
terms breaking the discrete $\mathbb{Z}_{2}$ symmetry between the SM and twin
sectors.

Let us first discuss scenario in which the only source of
$\mathbb{Z}_{2}$-breaking are the scalar couplings $\sigma$ and $\rho$. In
such a case, one may compute the tree-level $\mathbb{Z}_{2}$-asymmetric
contribution to the thermal masses of the Higgs bosons by evaluating the
second order derivative of the thermal potential $\eqref{eq:V_therm}$ at the
origin. This way, one obtains $\zeta_{A}\approx\zeta_{B}+\rho/2$, which,
substituted to equation (14), yields

$\frac{T_{A}}{T_{B}}\approx\sqrt{(1-\frac{\sigma}{\lambda})(1-\frac{\rho}{2\,\zeta_{A}})}.$
(17)

Thus, initial direction of symmetry breaking depends on two ratios:
$\sigma/\lambda$ and $\rho/(2\zeta_{A})$. $\rho/(2\zeta_{A})$ is usually
smaller since the tree-level thermal correction coefficient to the SM Higgs
mass $\zeta_{A}$ is greater than $\lambda$. Hence, it is mainly the sign of
$\sigma$ that determines whether $T_{A}$ or $T_{B}$ is larger (through the
first factor on the rhs of eq. (17)) so for $\sigma<0$, $h_{A}$ first acquires
non-zero VEV 666By definition, $\zeta_{A}$ is the coefficient of $T^{2}$ in
the SM Higgs thermal mass correction $\Pi_{h_{A}}(T)$. Looking at the
expression for $\Pi_{h_{A}}(T)$ given in the appendix B.2 one can verify that
for major part of considered parameter space $\lambda<\zeta_{A}$. Still, for a
small range of negative $\sigma$, for which
$|\rho|/(2\zeta_{A})\gtrsim|\sigma|/\lambda$, $h_{B}$ first acquires non-zero
VEV..

The twin symmetry between the model sectors can be also broken by couplings
which do not enter the tree-level scalar potential (2). Any significant
$\mathbb{Z}_{2}$-breaking interactions would generate hierarchy between
$\zeta_{A}$ and $\zeta_{B}$, which could overpower the effect of
$\sqrt{1-\sigma/\lambda}$ factor in eq. (14). In fact, by adjusting
$\zeta_{A}$ to $\zeta_{B}$ ratio, one gains control over the phase evolution
onset, irrespectively of $\mathbb{Z}_{2}$-breaking parameter values in the
scalar potential. On the other hand, strong $\mathbb{Z}_{2}$ breaking in the
non-scalar sector would generically induce large radiative corrections to the
SM Higgs mass parameter. This will obviously reintroduce fine-tuning to the
model unless one considers a UV extension where those radiative corrections
are naturally cancelled (an example of such extension is introduced in section
6).

In the following sections we describe the results of our numerical analysis.
We classify which paths of the thermal evolution of the fields may be present
depending on the parameters with a special emphasis of the paths in which
strong first-order phase transitions are present.

##  4 $\mathbb{Z}_{2}$ breaking in the Higgs potential only

Benchmark point | $f_{0}$ | $\lambda$ | $\kappa$ | $\sigma_{\text{tree}}+(\sigma_{\text{loop}})$ | $\rho_{\text{tree}}+(\rho_{\text{loop}})$
---|---|---|---|---|---
P1 | $4\,v_{A}$ | $1$ | $-0.1$ | $-0.12+(0.03)$ | $0.35+(-0.07)$
P2 | $4\,v_{A}$ | $1$ | $0.1$ | $0.083+(0.024)$ | $-0.054+(-0.072)$
Table 1: Benchmark points in TH model with $\mathbb{Z}_{2}$ symmetry broken
only in the scalar potential. The loop corrections to the tree level $\sigma$
and $\rho$ values are given in brackets.

###  4.1 Paths of thermal evolution

We start with a scenario in which $\mathbb{Z}_{2}$ breaking is obtained only
via the soft and hard breaking parameters in the tree-level TH potential. In
this class of models, two types of paths for the field evolution can be
present, which we discuss using two representative benchmark points listed in
table 1. We have numerically traced the thermal evolution of a global minimum
of the one-loop resummed effective potential (5). The resulting paths in
$h_{A}-h_{B}$ space together with the information about all the phase
transitions is shown in figure 1. As expected from the analytic discussion in
the previous section, we obtained different trajectories for the points with
positive and negative $\sigma$.

For $\sigma>0$ we found that both the path shape and phase transition order
are consistent with predictions made on analytical grounds in ref.
Fujikura:2018duw . In this regime, the symmetry is initially broken in the
$h_{B}$ direction and scalar phase does not undergo any tunneling – all phase
transitions were continues, see figure 1(a).

On the other hand, for non-negligible negative values of $\sigma$ 777Note that
the benchmark point with $\sigma<0$ is characterised by $\rho>0$ and
$\kappa<0$ which is a consequence of the constraints on the parameter space
from the minimization conditions and the measured Higgs mass of 125 GeV. we
found a new type of paths, which was not yet described in the literature in
context of TH models and is shown in figure 1(b). Here, the symmetry breaking
starts in the $h_{A}$ direction and then, as the universe cools down, the
$h_{A}$ accumulates its expectation value while $h_{B}$ remains zero. We found
that at some point global minimum exchanges with the local one, separated from
it by the potential barrier. Such behavior strongly indicates tunneling from
the local to the global minimum. During this phase transition $h_{A}$ moves to
zero and the EW symmetry is restored while $h_{B}$ gets a VEV and breaks twin
EW symmetry. Eventually, due to high twin Higgs VEV, masses of strongly
coupled twin particles become larger than the temperature of the Universe and
decouple from the thermal bath. At this point one recovers smooth one-
dimensional dynamics with $h_{A}$ behaving as Higgs field in unmodified SM and
$h_{B}$ fixed by the effective-field-theory (EFT) approximation

$h_{B}\approx\sqrt{f_{0}^{2}-h_{A}^{2}}.$ (18)

The scalar phase evolution proceeds in this manner until both, $h_{A}$ and
$h_{B}$, roll down to their zero temperature VEVs.

The parameters characterizing both phase transitions are given in table 2. To
ensure that in the case with $\sigma<0$ the FOPT really takes place, we have
numerically found the full instanton solution for the profile of the bubble of
the new phase (thus providing a simple existence prove). The aforementioned
solution was obtained with the so-called path deformation algorithm
implemented in the `cosmoTransitions` package Wainwright:2011kj .

Finally, let us comment on the accuracy of our one-loop computations. At both
benchmark points the condition (11) was not satisfied in the vicinity of the
phase transition. Therefore, we expect that for chosen points higher loop
terms in the effective potential were non-negligible at $T\sim T_{c}$ and thus
the one-loop results are merely qualitative estimates of the ongoing
processes. As we will show in the next section, for different parameter values
one can obtain larger $h_{I}(T_{c})/T_{c}$ which renders our results more
reliable.

$\sigma$ | tr. type | $T_{c}\;\text{[GeV]}$ | $T_{n}\;\text{[GeV]}$ | $v_{A}^{\text{ini}}\;\text{[GeV]}$ | $v_{B}^{\text{ini}}\;\text{[GeV]}$ | $v_{A}^{\text{fin}}\;\text{[GeV]}$ | $v_{B}^{\text{fin}}\;\text{[GeV]}$
---|---|---|---|---|---|---|---
$0.11$ | continues | $183$ | – | $0$ | $946$ | $0$ | $946$
$-0.09$ |  $1$-st order | $1235$ | $1172$ | $0$ | $652$ | $704$ | $0$
Table 2: Parameters characterizing phase transition obtained for benchmark
points from table 1 (for visualization see diagrams in figure 1). For the case
with negative $\sigma$ we have computed the full instanton solution and used
it to obtain the temperature at which the bubbles of the new phase start
nucleating which we refer to as nucleation temperature $T_{n}$. As one can
see, it only slightly differs from the critical temperature.
$v_{I}^{\text{ini}}$ ($v_{I}^{\text{fin}}$) corresponds to the VEV of $h_{I}$
at the beginning (end) of the phase transition.

(a) Phase diagram for benchmark point P1.  (b) Phase diagram for benchmark
point P2.

Figure 1: Evolution of the global minimum of the effective potential (5) for
the two benchmark points in table 1 with the opposite signs of $\sigma$. On
the right panel dashed arrow connects two degenerate minima, separated by the
potential barrier.

###  4.2 Systematic scan over the parameter space

As we have shown in the previous section, at some benchmark points the scalar
fields would likely undergo tunneling between the distant potential minima.
Such a process could have significant phenomenological consequences. If the
FOPT turns out to be sufficiently strong, it may open a way for EW-like
baryogenesis and/or generation of stochastic gravitational waves background
within TH. We are now going to analyze the model parameter space, classifying
the scalar phase transition order and its strength both in the SM and twin
sectors of the theory. Hereafter, we employ the standard condition for the
strong transition which originates from the requirement that needs to be
satisfied to decouple SM EW sphalerons888One may wonder if the standard
condition for sphaleron decoupling derived for the SM-like models is also
applicable to TH. We have checked, that due to structure similarities between
SM and TH, this condition should be roughly the same for both theories.

$\frac{\Delta v_{I}(T_{n})}{T_{n}}>1,$ (19)

where $\Delta v_{I}$ denotes the difference between the Higgs $I$ VEVs at the
end and at the beginning of the tunneling while $T_{n}$ is the temperature at
which stable bubbles of the new phase nucleate. We found that in the context
of TH the temperature of bubble nucleation is always well approximated with
the critical temperature $T_{c}$ (see for example table 1). Since computation
of the latter is much faster, in all the scans we use $\Delta
v_{I}(T_{c})/T_{c}$ rather than $\Delta v_{I}(T_{n})/T_{n}$.

If one allows for the modifications of the standard cosmological scenario by
allowing for the new relativistic species active at high temperatures, the
Hubble expansion rate during phase transition could be significantly boosted.
This weakens the sphaleron decoupling condition, down to Beniwal:2017eik

$\frac{\Delta v_{I}(T_{c})}{T_{c}}>0.5.$ (20)

The ratio $\Delta v_{I}(T_{c})/T_{c}$ used to determine the EW sphalerons
activity at phase transition encodes also an information about the quality of
the perturbative approximation (see the condition (11)). The larger is this
ratio in both sectors, the more robust are results obtained with the one-loop
effective potential. Thus, $\Delta v_{I}(T_{c})/T_{c}$ quantifies both the
strength and reliability of the phase transition and is the natural parameter
to be computed in the scans over the model parameter space.

It turns out that the condition for sphaleron decoupling in modified
cosmological scenarios (20) is slightly weaker than the perturbativity
condition (11). Nevertheless, rough character of the latter leaves the
possibility that the results with $0.5<\Delta v_{I}(T_{c})/T_{c}\lesssim g$
may still lead to approximately correct results. Since we would like to be as
general as it is possible, we are going to confront our results with both the
ordinary and weak sphaleron decoupling conditions.

In order to check where strong and weak conditions for FOPT are met, we run
the scan over the whole parameter space. Our program traced global minimum of
the effective potential (5) for approximately $2\,500$ points in the parameter
space. The minimum tracing was handled with the help of dedicated numerical
package `CosmoTransitions` which provides an accurate and fast solution for
this task Wainwright:2011kj .

(a) Transition strength in the SM sector (b) Transition strength in the TS
sector

Figure 2: The strength of the phase transition in the TH model with
$\mathbb{Z}_{2}$ broken only in the Higgs potential in $\kappa$-$\lambda$
plane for $f_{0}=4\,v_{A}$. The results for the SM (twin) Higgs are presented
on the left (right) panel. The regions where the SM (twin) Higgs transition
satisfies (19) is marked with dark blue (green) while the regions where
transition fulfills (20) is plotted with cyan (lime). In the red region the
transition was found very weak first-order or continues. The yellow region is
excluded because the $125\;\text{GeV}$ scalar corresponds to twin instead of
SM Higgs. The region above the green, blue and orange curves satisfy the
corresponding bounds $|\kappa|<\lambda$, $|\rho|<\lambda$, and
$|\sigma|<\lambda$. Black, vartical, dashed line corresponds to $\rho=0$.

The results of our computations are plotted in figure 2. It follows from the
plots that there exists a small region in the parameter space, where ordinary
sphaleron condition (19) for twin baryogenesis is satisfied 999Note that
sphaleron decoupling conditions for the SM EW baryogenesis are never
fulfilled. The ratio $v(T_{c})/T_{c}$ was computed using initial value of the
SM Higgs VEV since the final value is obviously close to zero. The ratio
obtained this way can only be used to roughly estimate the quality of
perturbative approximation.. This region is described by $\kappa<0$ which
follows from the fact that first-order phase transition in this scenario
requires $\sigma<0$ and hence relatively large positive $\rho$ which, in turn,
requires $\kappa<0$ to accommodate the 125 GeV Higgs mass. Even though non-
equilibrium dynamics takes place in both, the twin and SM sectors, in the
latter EW sphalerons do not freeze-out in the new phase and therefore the
formation of baryon asymmetry is not possible. Nevertheless, the obtained
results prove that phase transition in TH models can be strongly first-order.

To ensure that observed strong transitions are not an artifact of our low-
order perturbative expansion, we recomputed the ratio $v(T_{c})/T_{c}$, in the
model where all massive vectors are absent. Since gauge bosons are mainly
responsible for bad convergence at $T\sim T_{c}$, removing them renders theory
far more stable at phase transition. Even after removing gauge boson
contributions, which generically boost the transition order, at considerable
part of the parameter space FOPT remains sufficiently strong to fulfill the
condition (19). This proves that, at least for some parameter values, FOPTs
observed in TH do not require potentially unreliable gauge boson contributions
and are therefore robust.

We should, however, emphasize that the region where FOPT is strong is
characterized by relatively large values of $\rho$. For example, for
$\lambda=1$ we found that phase transition in the twin sector fulfills the
ordinary sphaleron decoupling condition (19) for $\rho\gtrsim 0.55$. Even
though this satisfies $\rho<\lambda$ so the 125 GeV can be described as PNGB
to a good approximation, some amount of tuning between $\kappa$ and $\rho$ is
required to obtain the measured Higgs mass for such large values of $\rho$.
This is because the 125 GeV Higgs mass implies $2\kappa+\rho\approx 0.1$ so
values of $\rho\gg 0.1$ require additional tuning. One cannot precisely
quantify this tuning without specifying UV completion of the model but one
should keep in mind that generically there will be some upper bound $\rho$
stemming from naturalness. The weaker twin sphaleron condition (20) for
$\lambda=1$ leads to $\rho\gtrsim 0.27$ for which the tuning in the Higgs mass
is less severe.

In conclusion, our results indicate that strong first-order transition, as
defined in (19), can take place in the ordinary TH model but may require some
tuning of the Higgs mass, whereas the weaker FOPT satisfying (20) can be
combined with substantial reduction of the tuning of the EW scale by the TH
mechanism.

##  5 Explicit $\mathbb{Z}_{2}$ breaking in lepton Yukawa couplings

In the previous section we assumed exact mirror $\mathbb{Z}_{2}$ symmetry
between the SM and TS couplings, broken only by the scalar interactions.
However, to alleviate the little hierarchy problem, it is sufficient to make
TS gauge couplings and the top Yukawa roughly equal to their SM analogues. To
avoid large radiative corrections to the twin Higgs mass parameter that would
result in fine-tuning one should also impose $\hat{y}_{q}\ll y_{t}$ for
$q\in\\{\hat{u},\hat{d},\hat{s},\hat{c},\hat{b}\\}$ (here, we use the hat to
denote TS fermions and their couplings). At first glance, it may seem that to
preserve the naturalness twin lepton Yukawas should also obey $\hat{y}_{l}\ll
y_{t}$ Craig:2015pha . However, if the cut-off scale for the lepton sector is
sufficiently low, which may still be allowed by the LHC searches particle2022
for lepton partners, large twin lepton Yukawa couplings may not lead to
conflict with the naturalness. Thus, in the present section we investigate the
effect of large twin Yukawa couplings on phase transitions in TH models. In
section 6 we will consider the case of large twin lepton Yukawa couplings in
the framework of SUSY TH model where it does not lead to fine-tuning provided
that sleptons are light enough Badziak:2019zys .

The possibility of the mirror symmetry breaking is often exploited in the
model variants that aim to resolve the problem of excessive dark radiation
which generically appears in TH models, see e.g. Chacko:2016hvu . The most
explicit example is the Fraternal TH in which TS contains only third
generation fermions Craig:2015pha . Another solution to the dark radiation
problem is to increase twin lepton and quark Yukawa couplings in order to
shift the dark QCD phase transition above the TS decoupling temperature
Barbieri:2016zxn ; Barbieri:2017opf .

To get a better understanding of how the discrete twin symmetry breaking in
lepton Yukawa(s) could affect the type and strength of phase transition, one
should look again at equation (14). Let us recall that in TH models strong
FOPT is feasible only when $T_{A}>T_{B}$. Previously, we considered the
specific case with $\zeta_{A}\sim\zeta_{B}+\rho/2$ and thus, strong FOPT
required $\sigma<0$. Negative $\sigma$ is always obtained at the cost of
introducing non-negligible hard $\mathbb{Z}_{2}$ breaking $\rho$. As it is
usually problematic to find UV completion of the TH model which naturally
explains big values of $\rho$ (see however Katz:2016wtw ), here we consider an
alternative approach that makes FOPT possible. It is based on manipulating
twin lepton Yukawa couplings in a way that introduces hierarchy between
$\zeta_{A}$ and $\zeta_{B}$. Since the SM Yukawa couplings for leptons are
negligible, the question really is whether introducing twin leptons with large
Yukawa couplings would boost the coefficient $\zeta_{B}$.

To address this question one can compute leptonic contributions to the twin
Higgs thermal mass $\zeta_{B}T^{2}$. In general, fermions can affect
$\zeta_{B}$ only through the thermal part of the effective potential (7). At
high temperatures, one can expand fermion thermal function
$J_{\text{F}}(\hat{m}_{l}^{2}/T^{2})$ (see Appendix B.1), to find the
correction from lepton species $i$ to the Higgs mass

$\Delta\hat{m}^{2}_{i\rightarrow h_{B}}=\frac{\partial^{2}}{\partial
h_{B}^{2}}\left[-\frac{\hat{n}_{i}T^{4}}{2\pi^{2}}J_{F}(\frac{\hat{m}_{i}^{2}}{T^{2}})\right]=T^{2}\frac{\hat{n}_{i}\hat{y}^{2}_{i}}{48}+\mathcal{O}\Big{(}\big{(}\frac{\hat{m}^{2}_{i}}{T^{2}}\big{)}^{\frac{3}{2}}\Big{)}$
(21)

where $\hat{n}_{i}$ is the number of degrees of freedom coupled with
$\hat{y}^{2}_{i}$. As one can see, twin leptons provide positive contribution
to $\zeta_{B}$, proportional to their Yukawa couplings squared and the number
of degrees of freedom. Whenever lepton couplings (or their number) are
sufficiently large, the symmetry is broken in the $h_{A}$ direction first, so
the necessary condition for FOPT can be satisfied even without introducing
$\rho\neq 0$.

At lower temperatures, the impact of twin fermions on TH phase transitions can
be understood with the help of EFT approximation (18). As long as high-
temperature expansion of the thermal function $J_{\text{F}}$ holds, the
leading order fermion corrections to the effective potential are of the form

$\displaystyle\delta V_{\text{eff}}=\left.\delta
V_{\text{eff}}\right|_{T=0}+\frac{1}{24}\sum_{i\in\text{Fermions}}n_{i}(m^{2}_{i}+\hat{m}^{2}_{i})+\mathcal{O}\Big{(}\big{(}\frac{m^{2}_{i}}{T^{2}}\big{)}^{\frac{3}{2}}\Big{)}=$
(22) $\displaystyle\left.\delta
V_{\text{eff}}\right|_{T=0}+\frac{1}{48}\sum_{i\in\text{Fermions}}\big{[}n_{i}(y^{2}_{i}-\hat{y}^{2}_{i})h_{A}^{2}+f_{0}^{2}\hat{y}^{2}_{i}n_{i}\big{]}+\mathcal{O}\Big{(}\big{(}\frac{m^{2}_{i}}{T^{2}}\big{)}^{\frac{3}{2}}\Big{)}.$

Thus, twin fermions with large Yukawa couplings will oppose EW symmetry
restoration. Usually, their effect becomes dominant at temperatures larger
than temperature of EW symmetry restoration $\sim 100\;\text{GeV}$, and hence,
they break it again pushing global minimum from $h_{A}=0$, $h_{B}=f_{0}$
towards the phase with $h_{A}=f_{0}$ and $h_{B}=0$.

The impact of heavy fermions on the phase transition in the models where the
SM Higgs is PNGB was noticed in Matsedonskyi:2020mlz . It was shown that for
sufficiently large number of strongly coupled twin fermions the EW symmetry is
not restored and the sphaleron decoupling condition (19) holds until $T\sim
f_{0}$. However, the number of fermions necessary to get symmetry non-
restoration together with sphaleron freeze-out exceeds the number of twin
leptons, even after including twin neutrinos101010To get the EW symmetry non-
restoration within TH it is necessary to introduce twin leptons and quarks
with order one Yukawas. Nonetheless, this is possible but requires very
specific UV completion to make it consistent with experimental bounds and keep
the fine-tuning under control Matsedonskyi:2020kuy .. Hence, in this section
we do not expect symmetry non-restoration.

Let us consider two cases. First: with the twin tau Yukawa coupling strongly
enhanced and negligible Yukawa couplings of lighter twin leptons. Second: with
all twin lepton Yukawas being large, which makes them significantly heavier
than their SM counterparts. We assume, for simplicity, that all twin lepton
Yukawa couplings are equal in this case.

(a) The SM transition strength.  (b) The TS transition strength.

Figure 3: The phase transition strength $\Delta v_{I}(T_{c})/T_{c}$ as a
function of twin lepton Yukawa coupling $\hat{y}_{l}$ for $\lambda=1$ and
$\rho=0$. Blue (orange) lines correspond to the case with the Yukawa couplings
of one (all three) twin leptons equal to $\hat{y}_{l}$. Dashed and solid lines
correspond to $f_{0}=3v_{A}$ and $f_{0}=4v_{A}$, respectively.

Let us start with a discussion of the phase transition strength dependence on
twin lepton Yukawa. We focus on the benchmark points with $\lambda=1$ and
$\rho=0$ and two different values of the $SU(4)$ breaking scales:
$f_{0}=3\,v_{A}$ and $f_{0}=4\,v_{A}$. The phase transition path is similar to
the one presented in figure 1(b). In figure 3 we present the plot of $\Delta
v_{I}(T_{c})/{T_{c}}$ against the enhanced twin lepton Yukawa coupling. We see
that the weak twin sphaleron decoupling condition $v_{B}(T_{c})/T_{c}>0.5$ is
satisfied for $\hat{y}_{l}$ at least about 0.4 (0.7) if all three (one) twin
lepton Yukawa couplings are large. We also see that in each case there is some
maximal value of $\hat{y}_{l}$ above which the strength of the phase
transitions decreases with $\hat{y}_{l}$. Whether the standard twin sphaleron
decoupling condition $v_{B}(T_{c})/T_{c}>1$ is satisfied depends on $f_{0}$.
This condition is satisfied for $f_{0}=4v_{A}$ and $\hat{y}_{l}$ at least
about 0.7 (1.1) if all three (one) twin lepton Yukawa coupling is large. The
sharp drop of the strength of the phase transition at some value of
$\hat{y}_{l}$ is due to an appearance of an intermediate phase so the fields
do not tunnel directly from a minimum with $(h_{A}\neq 0,h_{B}=0)$ to the one
with $(h_{A}=0,h_{B}\neq 0)$. This drop of the phase transition strength
appears for smaller values of $\hat{y}_{l}$ if $f_{0}$ is smaller, as clearly
seen in figure 3.

To visualize the combined effect of $\mathbb{Z}_{2}$-breaking in the TH scalar
potential and asymmetric Yukawa couplings we performed a scan of parameters
$\lambda$ and $\kappa$ which includes also the case of $\rho\neq 0$. In the
scan we considered the case in which only one twin lepton Yukawa coupling is
non-negligible and fixed to $\hat{y}_{l}=1.34$, which corresponds to the
maximal strength of the phase transition for $f_{0}=4v_{A}$, according to
figure 3. It is clearly seen in figure 4 that in such a case much larger parts
of parameters lead to FOPT, as compared to the case with no $\mathbb{Z}_{2}$
breaking in Yukawa couplings, cf. figure 2. Moreover, the region of parameters
in which the strong twin sphaleron decoupling condition is satisfied is now
characterized by positive values of $\kappa$ while $\rho$ can have both signs
or vanish. For $\rho=0$ strong FOPT corresponds to values of $\lambda$ in the
range between about $0.6$ and $1.5$ but one should keep in mind that this
range depends on the value of the twin lepton Yukawa couplings.

Figure 4: The same as in figure 2 but for the case of $\mathbb{Z}_{2}$
breaking by one twin Yukawa coupling $\hat{y}_{\tau}=1.34$ corresponding to
the maximal strength of the phase transition for the benchmark point with
$\lambda=1$, $\rho=0$ and $f_{0}=4v_{A}$, as demonstrated in figure 3.

##  6 Twin Higgs and lepton superpartners

The Twin Higgs model is not UV complete, and in principle should be treated as
an effective field theory limit of some high-energetic theory which solves the
big hierarchy problem by introducing colored top partners. In the present
work, we focus on SUSY TH models. We do not specify the mechanism that
generates the SU(4) invariant quartic term $\lambda$ which can originate
either from an $F$-term of a heavy singlet Falkowski:2006qq ; Chang:2006ra ;
Craig:2013fga  or a $D$-term of a new gauge symmetry badziak2017S ;
badziak2017M ; badziak2017A  and keep $\lambda$ as a free parameter. In this
effective description, the particle content of the model is that of MSSM and
twin MSSM.

If the superpartners are light enough, their existence can affect the scalar
field thermodynamics at TH energy scales. In particular, the order and
strength of phase transitions could be altered. Solution to the hierarchy
problem in SUSY TH models generically require the soft mass of the MSSM
sparticles and their corresponding twins to preserve the mirror
$\mathbb{Z}_{2}$ symmetry. The LHC set strong lower bounds on masses of
colored MSSM states, of about 1 to 2 TeV (with the exact limit depending on
the MSSM mass spectrum) particle2022 . On the other hand, the LHC limits on
slepton masses are much weaker and their masses $\mathcal{O}(100)$ GeV can be
still consistent with the experimental data. For such light sleptons even
$\mathcal{O}(1)$ twin lepton Yukawa couplings, which allow for strong FOPT as
shown in the previous section, do not lead to excessive fine-tuning of the EW
scale.

###  6.1 Setup

In this section, we investigate whether light (twin) sleptons can have non-
negligible impact on phase transitions in TH models. We perform this analysis
in a simplified model in which the rest of SUSY particles are much heavier and
decouple before the onset of phase transitions. We also assume that MSSM-like
Higgs bosons and their twin counterparts are decoupled, so during phase
transitions only the SM-like and twin Higgs bosons are dynamical. In terms of
the MSSM Higgs fields, they are defined as

$h_{I}=H^{u}_{I}\sin\beta_{I}+H^{d}_{I}\cos\beta_{I}\,,$ (23)

where up-type Higgs doublet $H^{u}_{I}$ couples to up-type fermions while
down-type Higgs doublet $H^{d}_{I}$ couples to down-type fermions and
$\beta_{I}\equiv\arctan\left(\langle H^{u}_{I}\rangle/\langle
H^{d}_{I}\rangle\right)$.

Since the Yukawa couplings for the SM leptons are very small, ordinary
sleptons do not affect the phase dynamics. However, twin sleptons may have
non-negligible impact if twin lepton Yukawa couplings are large enough because
they modify the effective potential (5) via Coleman-Weinberg and thermal
corrections. Here, we provide the formulae for slepton masses in order to
qualitatively analyze their effect. The tree-level diagonal entries of their
mass matrix read

$\displaystyle\hat{m}^{2}_{\text{sl}\;R}=\hat{\mu}^{2}_{R}+\frac{1}{2}\tilde{y}_{l}^{2}h_{B}^{2}\cos^{2}{\beta_{B}}-\frac{1}{4}g^{\prime
2}h_{B}^{2}\cos(2\beta_{B}),$ (24)
$\displaystyle\hat{m}^{2}_{\text{sl}\;L}=\hat{\mu}^{2}_{L}+\frac{1}{2}\tilde{y}_{l}^{2}h_{B}^{2}\cos^{2}{\beta_{B}}-\frac{1}{8}(g^{2}-\,g^{\prime
2})h_{B}^{2}\cos(2\beta_{B}),$ (25)

where $\hat{\mu}^{2}_{L(R)}$ is the soft SUSY breaking mass of the left-handed
(right-handed) slepton twin,
$\tilde{y}_{l}^{2}\equiv\hat{y}_{l}^{2}(\tan^{2}\beta+1)$ is the
supersymmetric twin lepton Yukawa. Mass terms proportional to gauge couplings
stem from the D-terms in the SUSY Lagrangian. They were derived from the
general expression Martin:1997ns  , valid for all squarks and sleptons

$m_{\phi\text{
D-term}}^{2}=\frac{1}{4}\left(T_{3\phi}-Q_{\phi}\sin^{2}\theta_{W}\right)(g^{2}+g^{\prime
2})v_{I}^{2}\cos(2\beta_{I}).$ (26)

Here, $T_{3\phi}$ stands for the third component of the weak isospin while
$Q_{\phi}$ corresponds to the electric charge of the left-handed
supermultiplet to which the field $\phi$ belongs.

In general, one can also consider the non-zero mixing between sleptons. In
such a case, the most universal, thermally corrected mass matrix reads

$M=\begin{pmatrix}\hat{m}^{2}_{\text{sl}\;R}+\hat{\Pi}_{\text{sl}\;R}&\;\;\frac{1}{\sqrt{2}}\tilde{y}_{l}h_{B}M_{\text{mix}}\\\\[5.0pt]
\frac{1}{\sqrt{2}}\tilde{y}_{l}h_{B}M_{\text{mix}}&\;\;\hat{m}^{2}_{\text{sl}\;L}+\hat{\Pi}_{\text{sl}\;L}\end{pmatrix}.$
(27)

Here, $M_{\text{mix}}$ is the mixing parameter which is the combination of
MSSM free parameters
$M_{\text{Mix}}=A_{l}\cos(\beta_{I})/\tilde{y}_{l}-\mu\sin(\beta)$, $A_{l}$ is
the soft trilinear lepton coupling while $\mu$ is the supersymmetric Higgsino
mass parameter. $\hat{\Pi}_{\text{sl}}(T)$ denotes thermal corrections which
have to be included in order to account for the leading order infrared
divergences. We neglect, for simplicity, the mixing between sleptons which is
justified unless $\mu$ or $\tan\beta$ are large.

###  6.2 Phenomenological impact of sleptons

The general analysis of the scalar dynamics in TH augmented by light twin
sleptons is a rather complex task. In principle, there are fourteen new free
parameters (six soft slepton masses, three mixing parameters, three twin
lepton Yukawas, and two SUSY $\beta$ angles) whose values could affect the
evolution of Higgs phases. Luckily, an overview of possible slepton effects
does not require scanning over all of them. The dynamics of Higgs VEVs depends
in the similar way on large groups of the new parameters, so there is no point
in considering them separately. Furthermore, it appears that some parameters
have negligible overall effects on the considered processes.

The motivation for introducing light twin sleptons was to obtain strong FOPT
in the most promising region of TH parameter space, where fine-tuning is
minimal. Hence, till the end of this section we will stick to the benchmark
point $\lambda=1$, $\rho=0$ and $f_{0}=4\,v_{A}$. As one could see from the
previous scans over TH parameter space, the region where FOPT is strong is
usually wide and does not contain any holes. Thus, if strong FOPT would take
place in the natural region of the model, there is virtually no chance that it
will be overlooked.

While there are no experimental constraints on the absolute values of
$\tan\beta_{\text{SM}}$ and $\tan\beta_{\text{TS}}$ the cancellation of the
loop corrections to SM Higgs mass requires them to be roughly the same. For
simplicity, throughout the phenomenological analysis we set them equal. It is
important to admit that for low $\tan\beta$ values the Landau pole of the top
Yukawa (and its twin) would appear at relatively low energy scales. With the
help of full one-loop MSSM renormalization group equations (RGEs)
Martin:1997ns  we have numerically checked that the top Yukawa Landau pole is
above $10^{4}\;(10^{16})\;\text{GeV}$ for $\tan\beta\gtrsim 0.4\;(1.5)$. If
the top quark is charged under some extra gauge symmetry, as it is the case
e.g. in SUSY D-term TH models badziak2017S ; badziak2017M ; badziak2017A , the
top Yukawa coupling may be asymptotically free even for smaller values of
$\tan\beta$.

In figure 5 we demonstrate how sleptons affect scalar phase transition for
several representative examples. Here, we set $\tan\beta=2$ and assume
$\mathbb{Z}_{2}$ symmetry between the soft masses in the model sectors. We
revisit these assumptions later in this section.

We start with the simple TH variant where twin tau Yukawa is enhanced and the
right-handed stau with its twin are the only SUSY partners present at thermal
equilibrium. The sleptons soft masses are set to $90\;\text{GeV}$ to keep the
physical mass of stau above the LEP bound particle2022 ; ALEPH:2001oot ;
DELPHI:2003uqw ; OPAL:2003nhx ; L3:2003fyi . The LHC limits on the right-
handed stau are still typically weaker than those from LEP CMS:2022rqk .

We also consider the case where all twin lepton Yukawas are large and all
right-handed sleptons are in thermal equilibrium. For simplicity, we set all
the twin lepton Yukawas equal. This assumption should not significantly affect
phase dynamics, as all twin sleptons contribute to the effective potential in
exactly the same way. The LEP lower bounds for all slepton masses are very
close to those for stau particle2022 ; ALEPH:2001oot ; DELPHI:2003uqw ;
OPAL:2003nhx ; L3:2003fyi  so all soft masses were set to $90\;\text{GeV}$.
The LHC limits on the masses of the right-handed smuon and selectron are
stronger than those for the right-handed stau but can be easily evaded if the
mass splitting between sleptons and the lightest supersymmetric particle (LSP)
is smaller than about 50 GeV ATLAS:2019lff .

(a) The SM transition strength.  (b) The TS transition strength.

Figure 5: The phase transition strength $\Delta v_{I}(T_{c})/T_{c}$ versus
twin lepton Yukawa coupling $\tilde{y}_{l}$ in SUSY TH model with light
sleptons for $\tan\beta=2$. The blue and cyan lines correspond to the scenario
with enhanced $\tilde{y}_{l}$ for only one family of twin leptons, while the
pink and red lines correspond to Yukawa couplings of all three twin leptons
set to $\tilde{y}_{l}$. The cyan and pink lines denote the results for right-
handed sleptons only, wheres the red and blue lines correspond to the scenario
with both slepton chiralities. Solid lines denote the results with
$\hat{\mu}_{\text{R (L)}}=90\;\text{GeV}$ while dashed lines for those with
$\hat{\mu}_{\text{R (L)}}=290\;\text{GeV}$.

We analyze scenarios with both sleptons chiralities (but no mixing). The LHC
limits for degenerate left-handed and right-handed sleptons are much stronger
than the corresponding limits on the right-handed sleptons so in this scenario
sleptons with masses of about 100 GeV are allowed only for fine-tuned values
of the mass splitting between sleptons and the LSP ATLAS:2019lng ;
ATLAS:2022fxl . Thus, we present results not only for
$\hat{\mu}_{\text{sl}}\equiv\hat{\mu}_{\,\text{L}}=\hat{\mu}_{\text{R}}=90\;\text{GeV}$
but also for $\hat{\mu}_{\text{sl}}=290\;\text{GeV}$. In the latter case, the
physical masses of sleptons are around $300\;\text{GeV}$ so the LHC bounds are
satisfied as long as the mass splitting between sleptons and the LSP is below
about 100 GeV ATLAS:2019lff .

(a) The SM transition strength.  (b) The TS transition strength.

Figure 6: Contour plots of the phase transition strength $\Delta
v_{I}(T_{c})/T_{c}$ in the $\tilde{y}_{l}$-$\tan\beta$ plane in SUSY TH model
with light right- and left-handed sleptons from all three families. The soft
slepton masses were set to
$\mu_{\text{sl}}=\hat{\mu}_{\text{sl}}=100\;\text{GeV}$. In accordance with
MSSM RGEs the Landau pole of the top Yukawa appears below
$10^{16}\;\text{GeV}$ if $\tan\beta\lesssim 1.5$ (region under the red line).
We cut the plot at $\tan\beta=0.4$ since for smaller $\tan\beta$ values the
Landau pole would appear below $10^{4}\;\text{GeV}$. In the yellow region a
Landau pole for the twin lepton Yukawa coupling appears below $10\,f_{0}$. In
the blue region the MSSM slepton mass squared was found negative at some point
of the thermal evolution. In the green region sphaleron wash-out of SM baryon
asymmetry occurs only at $T>T_{c}$.

It can be seen in figure 5, that strong FOPT satisfying the condition (19) can
be met in all scenarios for the moderate $\tilde{y}_{l}$ values. In fact, the
results of the scan presented in figure 5 resemble the results in section 5
(keep in mind that $\tilde{y}_{l}=\hat{y}_{l}/\cos\beta$) although this
similarity is not exact due to slepton effects, D-terms and thermal mass
explicit dependence on $\tan\beta$. As one could expect, the transition type
is not affected when twin Yukawas are small because in this regime twin
leptons and their superpartners contribution to the effective potential is
negligible. Only when $\tilde{y}_{l}$ exceeds some critical value (which
depends on the scenario) its influence on the transition type becomes
important. Then, for larger $\tilde{y}_{l}$ values a temporal drop in $\Delta
v_{I}(T_{c})/T_{c}$ occurs. It is caused by the new intermediate phase that
appears between the minima previously separated by the potential barrier.

The sleptonic impact on phase transition order and strength was observed in
all considered scenarios. For very large twin lepton Yukawa couplings
$\tilde{y}_{l}\gtrsim 2$ which are close to the boundary of validity of
perturbative calculations, we observe strong FOPT induced by light twin
sleptons. Avoiding a Landau pole for the twin lepton Yukawa coupling below
$10\,f_{0}$ sets an upper bound on $\tilde{y}_{l}\lesssim 3$. We admit, that
for the allowed values of $\tilde{y}_{l}$ only in the scenario with three
families of light sleptons with large twin Yukawa couplings the phase
transition induced by twin sleptons can be stronger than that induced by twin
leptons.

Let us now discuss a dependence of the strength of the phase transition on
$\tan\beta$. For large $\tan\beta$ the SM-like Higgs boson and its twin are
strongly dominated by the corresponding up-type Higgs fields $H^{u}_{I}$. On
the other hand, twin sleptons couple to the down-type Higgs field $H^{d}_{B}$.
Therefore, the impact of twin sleptons on the strength transition decreases
with $\tan\beta$. This can be seen in figure 6 where we plot contours of
$\Delta v_{I}/T_{c}$ in the $\tilde{y}_{l}$-$\tan\beta$ plane for a benchmark
scenario with all three families of the right- and left-handed twin sleptons
with the soft slepton masses fixed to $100\;\text{GeV}$ 111111We choose this
scenario as it offers the strongest transition in the TS and thus gives best
prospects for DS baryogenesis and emission of the detectable wave
gravitational background (see section 7 for details). For small $\tan\beta$
the phase transition can become very strong with $v_{B}/T_{c}\gtrsim 2$ for
$\tan\beta\lesssim 1.5$ and perturbative values of twin lepton Yukawa
couplings. For example, for $\tan\beta\approx 1$ slepton-induced FOPT with
$v_{B}/T_{c}>1$ requires $\tilde{y}_{l}\gtrsim 1$ (in addition to the one
induced by twin leptons which occurs for $\tilde{y}_{l}\approx 0.5$). These
coupling values correspond to a region in parameter space in which the top and
twin tau Yukawa couplings are perturbative up to scales many orders of
magnitude above the EW scale.

We observe that when $\tan\beta<1$ the D-term contributions in the expression
for slepton masses (24) become negative. In the SM sector lepton Yukawas are
tiny and negative D-terms could dominate over other terms including thermal
mass, provided $h_{A}(T)$ is big. This is indeed the case in the blue region
in figure 6, where the physical masses of SM sleptons were found negative at
some $T>0$, which indicates spontaneous breaking of electromagnetic symmetry.

Figure 7: Left plot: $h_{A}(T)$ (blue line) and $h_{B}(T)$ (red line) VEVs
during thermal evolution. Right plot: path in the $\langle
h_{A}\rangle-\langle h_{B}\rangle$ space with the part of curve where SM
sphalerons are suppressed (i.e. $v_{A}(T)/T>1$) marked in green. The plots
were completed at the test point $\tilde{y}_{l}=2.7$, $\tan\beta=0.8$,
$\mu_{\text{sl}}=\hat{\mu}_{\text{sl}}=100$ GeV placed in the green region in
figure 6. Baryon asymmetry wash-out does not occur below $T\approx
800\;\text{GeV}$.

We also found that for large (but still perturbative) values of twin lepton
Yukawa coupling and $\tan\beta\lesssim 1$ the thermal evolution of the fields
is qualitatively different. This region is depicted in green in figure 6 and
the transition path for a representative point is displayed in figure 7. In
contrast to transition paths for larger $\tan\beta$, after the SM-like Higgs
gets a VEV at temperature around $f_{0}$, the EW symmetry is never restored
for lower temperatures. As a consequence, SM sphalerons are active only at
temperatures significantly above the EW scale. For the benchmark point in
figure 7 the SM sphaleron decoupling occurs around 800 GeV, so there is no
wash-out of the SM baryon asymmetry below that temperature. Such cosmological
evolution may have important implications for models of baryogenesis.

Figure 8: The phase transition strength $\Delta v_{I}(T_{c})/T_{c}$ vs twin
slepton soft mass squared $\hat{\mu}^{2}_{\text{sl}}$ (equal among all the
three generations of the right- and left-handed twin sleptons) for
$\tan\beta=2$ and $\mu_{\text{sl}}=90$ GeV. Blue, orange and green lines
correspond to $\tilde{y}_{l}$ set to $2$, $2.5,$ and $3$, respectively. Each
line ends at the point where one of the physical twin slepton masses was found
negative.

Let us now discuss how the strength of the phase transition depends on the
value of the soft twin slepton masses. In figure 8 we present $\Delta
v_{I}(T_{c})/T_{c}$ as a function of $\hat{\mu}_{\text{sl}}^{2}$ being a
common soft mass for all three families of the right- and left-handed twin
sleptons. We fix $\tan\beta=2$ and show the results for three choices of twin
lepton Yukawa couplings (equal among the generations): $2$, $2.5$ and $3$. We
see that the phase transition strength $\Delta v_{B}(T_{c})/T_{c}$ is anti-
correlated with $\hat{\mu}_{\text{sl}}^{2}$, leading to stronger transitions
for smaller soft twin slepton mass squared. On the other hand, when
$\hat{\mu}_{\text{sl}}^{2}\gtrsim f_{0}^{2}$ twin sleptons effectively
decouple, and we recover results obtained for non-supersymmetric TH with
enhanced twin lepton Yukawas.

In figure 8 we also see that $v_{B}(T_{c})/T_{c}$ grows significantly when
$\hat{\mu}_{\text{sl}}^{2}$ becomes negative. For large range of negative
values of $\hat{\mu}_{\text{sl}}^{2}$ the physical masses squared of twin
sleptons remain positive due to large positive supersymmetric contributions
enhanced by large values of $v_{B}$ and $\tilde{y}_{l}$. However, large
negative values of $\mu_{\text{sl}}^{2}$ make the MSSM sleptons tachyonic.
Hence, enhancing the strength of twin phase transition by negative values of
$\hat{\mu}_{\text{sl}}^{2}$ requires $\mathbb{Z}_{2}$-breaking in the soft
slepton masses. Such $\mathbb{Z}_{2}$-breaking generically occurs for large
twin lepton Yukawa coupling due to negative one-loop corrections
$\sim\tilde{y}_{l}^{2}$ to twin slepton soft masses, while the corresponding
corrections to MSSM sleptons are negligible. Thus, the slepton soft masses
consist of two parts: a positive $\mathbb{Z}_{2}$-preserving contribution and
the $\mathbb{Z}_{2}$-breaking negative loop corrections which, for
sufficiently large $\tilde{y}_{l}$, may become dominant. In figure 8 we set
$\mu_{\text{sl}}=90$ GeV in order to satisfy the experimental bounds on
slepton masses but $\Delta v_{B}(T_{c})/T_{c}$ does not depend much on a
specific value of $\mu_{\text{sl}}$.

In summary, the sphaleron decoupling condition (19) in the twin sector is
satisfied in all considered scenarios, at least for some values of
$\tilde{y}_{l}$. Even though the SM sphalerons never decouple in a way that
could lead to the EW-like baryogenesis in the SM sector, there is a region of
parameter space where baryon asymmetry wash-out does not occur. This opens a
way for baryon asymmetry transfer from the twin sector and/or inheritance of
baryon asymmetry produced in the earlier stages of the evolution of the
Universe.

##  7 Stochastic gravitational wave background

The strong FOPTs are often followed by emission of the stochastic wave
gravitational background. Since there are regions of the model parameter space
where such transitions are feasible, one may numerically simulate
corresponding GW spectra. If, for some parameter choice, the obtained signal
falls in the range of any existing or planned detector, one will obtain
experimental predictions that could be verified in the near future.

The accurate estimation of GW spectrum emitted during FOPT is a complex
problem. Here, we use the standard approximation Caprini:2015zlo ;
Caprini:2019egz , which is based on four input parameters: the temperature of
bubble percolation $T_{p}$ at which the probability of finding a point still
in the false vacuum is 0.7, latent heat released during phase transition
$\alpha$, its inverse duration time in the Hubble units
$\beta_{\ast}/H_{\ast}$, and the terminal velocity of the new phase bubble
walls $v_{w}$. Further in this section we approximate $T_{p}$ with the
temperature of critical bubble nucleation $T_{n}$, which is accurate for
$\alpha\lesssim 0.2$ Wang:2020jrd . To compute the necessary parameters one
has to find the profile of the new phase bubbles that nucleate during FOPT.
This is a rather complex numerical problem which was handled with dedicated
Python package `CosmoTransitions` with custom modifications Wainwright:2011kj
; marcopriv . More details about the computation of the input parameters and
method used for estimation of GW signal are given in appendix C.

line style | $\mathbf{f_{0}}$ | $\boldsymbol{\lambda}$ | $\boldsymbol{\kappa}$ | $\mathbf{T_{n}}$ | $\boldsymbol{\alpha}$ | $\mathbf{\boldsymbol{\beta}_{\ast}/H_{\ast}}$ | $\boldsymbol{v_{w}}$
---|---|---|---|---|---|---|---
solid red | $4\,v_{A}$ | $1$ | $-0.1$ | $1172$ | $1.3\times 10^{-4}$ | $5.1\times 10^{3}$ | $0.39$
dashed red | $4\,v_{A}$ | $1$ | $-0.25$ | $821$ | $9.4\times 10^{-4}$ | $1.5\times 10^{4}$ | $0.33$
solid blue | $8\,v_{A}$ | $1$ | $-0.1$ | $2171$ | $1.9\times 10^{-4}$ | $2.4\times 10^{3}$ | $0.51$
dashed blue | $8\,v_{A}$ | $1$ | $-0.25$ | $1422$ | $1.6\times 10^{-3}$ | $4.8\times 10^{3}$ | $0.54$
Table 3: Test points for the UV agnostic TH. The corresponding GW spectra are
plotted in figure 9(a).

The ratio $\Delta v_{I}(T_{c})/T_{c}$ used throughout this paper to estimate
phase transition strength does not appear in any expression for the GW
spectrum. Nonetheless, detectable wave gravitational background cannot be
generated when transition is not first-order.The one-loop computations are
inaccurate when $\Delta v_{I}(T_{c})/T_{c}$ is too low (see the appendix B.2)
and thus, the transitions that were found to be weak first-order at one loop
often turn out to be continues when more accurate computations are performed.
For the above reason, in this section we restrict our attention to the
parameter values for which the condition (11) is fulfilled.

In order to accurately identify regions in the parameter space where
potentially detectable gravitational wave background is generated, one should
run full scan over the parameter values for which FOPT occurs. Here, we
refrained from performing such scan, but instead computed the GW signal for
carefully chosen test points, listed in tables 3 and 4. For the UV agnostic TH
we use the free parameter values for which the strong (in terms of $\Delta
v_{I}(T_{c})/T_{c}$) FOPT happens but no significant tuning is present121212We
confirmed numerically that for the UV agnostic TH latent heat of the
transition $\alpha$ and hence also the GW signal is positively correlated with
the ratio $\Delta v_{I}(T_{c})/T_{c}$.. In the case of SUSY extension, we
restrict our attention to benchmark points with $\lambda=1$, $\rho=0$,
$f_{0}=4\,v_{A}$ and consider different values of twin lepton Yukawas, slepton
soft masses, and $\tan\beta$ values for two distinct scenarios introduced in
section 6.2. The corresponding spectra together with present and future
detector sensitivities are shown in figure 9.

line style | $\mathbf{\hat{n}_{l}}$ | tan$\boldsymbol{\beta}$ |  $\mathbf{\hat{\boldsymbol{\mu}}^{2}_{\text{sl}}}$ [GeV2] |  $\mathbf{T_{n}}$ [GeV] | $\boldsymbol{\alpha}$ | $\mathbf{\boldsymbol{\beta}_{\ast}/H_{\ast}}$ | $\mathbf{v_{w}}$
---|---|---|---|---|---|---|---
solid green | $2$ | $2$ | $90^{2}$ | $586$ | $3.3\times 10^{-3}$ | $2.3\times 10^{4}$ | $0.21$
dashed green | $2$ | $2$ | $-300^{2}$ | $565$ | $4.3\times 10^{-3}$ | $1.1\times 10^{4}$ | $0.30$
solid blue | $2$ | $1$ | $90^{2}$ | $514$ | $7.6\times 10^{-3}$ | $1.7\times 10^{4}$ | $0.23$
dashed blue | $2$ | $1$ | $-375^{2}$ | $491$ | $1.0\times 10^{-2}$ | $7.9\times 10^{3}$ | $0.34$
solid light green | $12$ | $2$ | $90^{2}$ | $402$ | $1.9\times 10^{-2}$ | $1.5\times 10^{4}$ | $0.24$
dashed light green | $12$ | $2$ | $-310^{2}$ | $375$ | $2.8\times 10^{-2}$ | $5.4\times 10^{3}$ | $0.38$
solid light blue | $12$ | $1$ | $90^{2}$ | $336$ | $5.0\times 10^{-2}$ | $5.5\times 10^{3}$ | $0.40$
dashed light blue | $12$ | $1$ | $-230^{2}$ | $317$ | $6.2\times 10^{-2}$ | $3.2\times 10^{3}$ | $0.49$
Table 4: Test points for the SUSY TH model with light sleptons and large twin
lepton Yukawa couplings. We have considered two scenarios–one with a single
twin lepton with enhanced Yukawa and right-handed slepton
($\hat{n}_{\text{sl}}=2$) and another, where all Yukawas of charged lepton
twins are big and sleptons with both chiralities are active
($\hat{n}_{\text{sl}}=12$). We fix $\tilde{y}_{l}=3$ since we have numerically
verified that in all considered points the GW signal is strongest for the
largest $\tilde{y}_{l}$. The remaining parameters are fixed at $\lambda=1$,
$\rho=0$, $f_{0}=4\,v_{A}$. The corresponding GW spectra are plotted in figure
9(b).

In the non-supersymmetric TH model the GW signal is many orders of magnitude
below the sensitivity of future detectors. The GW signal increases somewhat
with the $SU(4)$ breaking scale $f_{0}$. This is because the GW signal grows
with the transition temperature which is proportional to $f_{0}$.
Nevertheless, even for $f_{0}=8v_{A}$ the GW signal is more than five orders
of magnitude smaller than the sensitivity of future detectors.

In the supersymmetric extension of TH introduced in section 6 the transition
could be caused either by twin leptons or their SUSY partners. In the twin
lepton-induced case, which corresponds to the first maximum on the plots in
figure 5, we found the GW signal to be very weak. Therefore, in our numerical
analysis we have focused on the slepton-induced transition, which occurs close
to the second maximum on the plots in figure 5. We found that the most
promising case for generation large GW signal is the one in which all three
twin lepton Yukawa couplings are large, twin sleptons from all three families
are as light as possible and $\tan\beta$ is small. Even though the
corresponding GW signal is much stronger than in non-supersymmetric case, its
peak frequency is still more than an order of magnitude off from the reach of
AEDGE Bertoldi:2021rqk  and Einstein Telescope Punturo:2010zz ; Hild:2010id .
The obtained power-low integrated spectra fall in the blind spot of currently
planned detectors and for a rather large range of frequencies between
$\mathcal{10^{-1}}$ and $\mathcal{10^{1}}$ Hz is roughly two orders of
magnitude below their sensitivity. Nevertheless, to exclude the detection
possibility of this scenario more accurate estimation of GW spectrum would be
necessary.

(a) GW spectra in UV agnostic TH model. (b) GW spectra in SUSY TH model.

Figure 9: GW spectra for the benchmark points presented in tables 3 and 4.
Shaded regions mark the integrated power-law sensitivity Thrane:2013oya  of
the existing and future GW detectors: LIGO collaboration2015advanced , LISA
Bartolo:2016ami ; Caprini:2019pxz , AION Badurina:2019hst , AEDGE
Bertoldi:2021rqk  and Einstein Telescope Punturo:2010zz ; Hild:2010id . The
LISA sensitivity was plotted based on data from the web tool PTPlot
Caprini:2019egz . Data about the rest of detector ranges courtesy of M.
Lewicki lewickipriv .

##  8 Conclusions and Outlook

We have shown that Twin Higgs model features strong first-order phase
transition provided that appropriate source of $\mathbb{Z}_{2}$ symmetry
breaking between twin and ordinary matter sectors is present. Thermal
evolution of the Universe follows an unusual path in the $h_{A}$–$h_{B}$
plane, as seen in figure 1(b), since it is the SM-like Higgs which acquires a
non-zero expectation value first at temperature $\sim f_{0}$. As the Universe
cools down, the FOPT occurs during which the SM-like and twin Higgs change
their field expectation values. The VEV of the SM-like Higgs decreases from a
large initial value to a vanishing one, which excludes the possibility of EW-
like baryogenesis. On the other hand, the VEV of the twin Higgs increases,
leading to twin sphaleron decoupling after the transition completes. This
opens a possibility for dark baryogenesis that may induce the observed
asymmetry of SM baryons if an appropriate portal between the twin and SM
sector is present.

Our results also show that it is crucial to analyze two-field thermal
evolution of the Higgs fields to get a complete picture of phase transitions
in TH models. The previous study Fujikura:2018duw , performed using one-field
approximation, did not find first-order phase transitions.

We have proposed two distinct scenarios of $\mathbb{Z}_{2}$ symmetry breaking
that leads to the FOPT described above. In the first one positive, big
asymmetric quartic coupling is present in the tree-level scalar potential.
This way one can indeed obtain strong FOPT within the original TH model
without specifying its UV completion. Strong FOPT in this scenario requires
hard $\mathbb{Z}_{2}$ breaking, which leads to some tuning against negative
$\mathbb{Z}_{2}$-preserving but $SU(4)$-breaking quartic coupling to
accommodate the 125 GeV Higgs mass. Such region of the UV agnostic TH
parameter space has been often left without proper analysis.

The second possibility is to include explicit $\mathbb{Z}_{2}$ breaking by
increasing Yukawa couplings of twin leptons. Interestingly, $\mathbb{Z}_{2}$
breaking in Yukawa couplings relaxes also the dark radiation problem in TH
models Barbieri:2016zxn ; Barbieri:2017opf  and makes twin neutralino a
natural candidate for thermally produced DM Badziak:2019zys . In this
scenario, strong FOPT can be induced by twin leptons without introducing hard
$\mathbb{Z}_{2}$ breaking in the scalar potential if twin lepton Yukawa
couplings are around one.

In order to avoid excessive tuning in the Higgs potential stemming from large
twin lepton Yukawa couplings the TH model must be UV completed at relatively
low scale. We have achieved it by considering SUSY UV completion with light
sleptons. Moreover, we found that light twin sleptons lead to even stronger
FOPT provided that $\tan\beta\lesssim 2$ and twin lepton Yukawa couplings are
above one.

It is worth mentioning that in the SUSY TH model with light sleptons
characterized by small $\tan\beta\lesssim 1$ and very large (but still
perturbative) twin Yukawa couplings the EW symmetry is broken for temperatures
up to above 1 TeV, see figure 7 for a representative thermal field evolution
in such a case. In this scenario, wash-out of the SM baryon asymmetry, which
might be created at some earlier stages of the evolution of the Universe, is
delayed until $T\sim\text{TeV}$.

An important part of this work was devoted to estimation of the spectrum of
stochastic gravitational wave background which could be emitted during the
FOPT. We found that in the UV agnostic TH the obtained GW spectra are nowhere
near the reach of the present or planned detectors. In contrast, in the
supersymmetric TH model with light sleptons the gravitational wave signal from
FOPTs is stronger by many orders of magnitude. In a small region of parameter
space we found that the GW signal for frequencies within reach of AEDGE
Bertoldi:2021rqk  and Einstein Telescope Punturo:2010zz ; Hild:2010id
detectors can be close but still about two orders of magnitude below their
expected sensitivity.

Our findings open several new avenues for research. Strong FOPTs during which
twin EW sphalerons decouple may allow for generation of twin baryon asymmetry.
It would be interesting to investigate scenarios in which such asymmetry is
transferred to the SM sector and leads to the observed abundance of the SM
baryons. It would be also interesting to investigate other variants of TH
models in which $\mathbb{Z}_{2}$ breaking leads to strong FOPTs. Finding
models in which detectable GWs and/or EW symmetry non-restoration up to
temperatures much above the EW scale are more generic phenomena would be
particularly exciting.

## Acknowledgments

The authors would like to thank Keisuke Harigaya, Marek Lewicki, Marco
Merchand, Mariano Quiros, Bogumiła Świeżewska and Mateusz Zych for useful
discussions and Kohei Fujikura for useful correspondence. This work was
partially supported by the National Science Centre, Poland, under research
grant no. 2020/38/E/ST2/00243.

##  Appendix A Tree-level Higgs masses and VEVs beyond PNGB approximation.

The exact expressions for the tree-level masses and VEVs of SM and twin
Higgses can be derived analytically Ahmed:2017psb . They read

$v_{A}^{2}=\lambda
f_{0}^{2}\frac{-\sigma+\kappa(1-\frac{\sigma}{\lambda})}{\lambda\rho+\kappa(2\lambda+\rho+\kappa)},\qquad
v_{B}^{2}=f^{2}-v_{A}^{2}\,,$ (28)

where

$f^{2}=f_{0}^{2}\frac{\lambda\rho+\kappa(2\lambda-\sigma)}{\lambda\rho+\kappa(2\lambda+\rho+\kappa))}\,,$
(29)

The expressions for the tree-level masses of the SM-like Higgs boson, $m_{h}$,
and its twin counterpart, $m_{h^{\prime}}$, are

$\displaystyle m_{h(h^{\prime})}^{2}=\rho
v_{A}^{2}+f^{2}(\lambda+\kappa)\left(1\mp\sqrt{1-A}\right)\,,$ (30)
$\displaystyle A\equiv
2\frac{v_{A}^{2}}{f^{2}}\frac{\lambda\rho+\kappa(4\lambda+\rho+2\kappa)}{(\lambda+\kappa)^{2}}-\frac{v_{A}^{4}}{f^{4}}\frac{4\lambda\rho+\rho^{2}+\kappa(8\lambda+4\rho+4\kappa)}{(\lambda+\kappa)^{2}}\,,$

where the upper (lower) sign in the above equation corresponds to $m_{h}$
($m_{h^{\prime}}$). The approximate expressions (3) and (4) are recovered by
expanding the exact tree-level formulas to the leading order in
$v_{A}^{2}/f^{2}\approx v_{A}^{2}/f^{2}_{0}$ and
$\kappa/\lambda,\;\sigma/\lambda,\;\rho/\lambda$.

##  Appendix B Finite temperature terms in the effective potential

###  B.1 Low and high temperature expansion of the thermal potential

To evaluate thermal functions $J_{B}$ and $J_{F}$ one can directly compute the
integrals used in their definitions (8). Alternatively, it is possible to use
the expressions for the low and high temperature expansion of those functions.
For $T\gg m$ one can use Quiros:1999jp

$\displaystyle
J_{B}(\frac{m^{2}}{T^{2}})\approx-\frac{\pi^{4}}{45}+\frac{\pi^{2}}{12}\frac{m^{2}}{T^{2}}-\frac{\pi}{6}\left(\frac{m^{2}}{T^{2}}\right)^{\frac{3}{2}}-\frac{m^{4}}{32T^{4}}\log\left(\frac{m^{2}}{a_{b}\,T^{2}}\right)+\mathcal{O}\left(\frac{m^{6}}{T^{6}}\right)$
(31) $\displaystyle
J_{f}(\frac{m^{2}}{T^{2}})\approx\frac{7\pi^{4}}{360}-\frac{\pi^{2}}{24}\frac{m^{2}}{T^{2}}-\frac{m^{4}}{32T^{4}}\log\left(\frac{m^{2}}{a_{f}\,T^{2}}\right)+\mathcal{O}\left(\frac{m^{6}}{T^{6}}\right)$
(32)

with $a_{b}=16\,\pi^{2}e^{\frac{3}{2}-2\gamma_{\text{E}}}$ and
$a_{f}=\pi^{2}e^{\frac{3}{2}-2\gamma_{\text{E}}}$, where $\gamma_{\text{E}}$
is the Euler’s constat.

For $T\leq m$ the corresponding expansion reads Curtin:2016urg

$\displaystyle
J_{B}(\frac{m^{2}}{T^{2}})=-\sum_{n=1}^{\infty}\frac{1}{n^{2}}\frac{m^{2}}{T^{2}}\text{K}_{2}(n\,\frac{m}{T})$
(33) $\displaystyle
J_{F}(\frac{m^{2}}{T^{2}})=-\sum_{n=1}^{\infty}\frac{(-1)^{n}}{n^{2}}\frac{m^{2}}{T^{2}}\text{K}_{2}(n\,\frac{m}{T})$

where $\text{K}_{2}(x)$ is the Bessel function of the second kind. The above
expression for low temperature expansion is less intuitive but it remains
surprisingly accurate even for relatively high temperatures. It also does not
require including many terms to get a good precision – two or three are
usually sufficient.

The computation of the approximate forms (31)-(33) is much faster than direct
evaluation of the integrals (8). Since in the whole range $J_{B}$ and $J_{F}$
are well-approximated either by low or high temperature expansion, to obtain
numerical results we use the series defined above rather than $J_{B}$ and
$J_{F}$ definition.

###  B.2 Resummation of infrared divergences

At high temperatures the non-abelian gauge theories with massive vectors
suffer from infrared divergences. The divergent contributions arise in the
vicinity of phase transition, where all particles which acquire mass due to
Higgs mechanism become massless. The leading order infrared contribution stems
from the so-called daisy diagrams. It scales with the number of massless
bosonic loops $n$ as $\sim(\mathcal{g}^{2}T^{2}/m^{2})^{n}$, where
$\mathcal{g}$ is the characteristic gauge coupling and $m$ characteristic
gauge mass scale of the problem. At the phase transition
$m\approx\mathcal{g}\phi\sim\mathcal{g}T$ and the scaling factor becomes
roughly $1$.

This issue can be alleviated by the resummation of the leading order infrared
divergences. After the resummation, the series of most divergent diagrams is
suppressed by an additional factor
$\mathcal{g}^{2}T/m(\phi)\sim\mathcal{g}v_{c}/T_{c}$131313In principle, in the
suppression factor one should use thermally resummed mass $\overline{m}(\phi)$
defined by the eq. (35). However, most divergent contributions come from
transverse modes of the massive vectors which do not receive thermal
corrections at leading order, and thus $\overline{m}(\phi)\approx m(\phi)$..
Hence, the rough criterion for the applicability of one-loop approximation of
the effective potential after the resummation reads

$\frac{\phi(T)}{T}>\mathcal{g},$ (34)

from which one obtains the condition (11) for TH.

Strictly speaking, the results derived within one-loop approximation in gauge
theories are reliable only when the transition is sufficiently strong, i.e.
$\phi(T_{c})/T_{c}\gg\mathcal{g}$ (otherwise, the computations are unreliable
in the vicinity of phase transition). Moreover, the accuracy increases with
the transition strength. To some extent, the criterion (34) may be weakened by
the more accurate resummation techniques Curtin:2016urg  which take into
account the sub-leading multi-loop corrections. However, their implementation
requires much more effort and usually does not improve the accuracy a lot at
the temperature at which the strong FOPT occurs.

In this paper, we consider two fields dynamics, and hence one field criterion
(11) does not necessarily guarantee the perturbativity at phase transition. Up
to our knowledge, the proper criterion for the type of FOPT we found in TH is
yet unknown. We expect that if the condition (11) is satisfied in the theory
sector with broken symmetry before and after the transition, our perturbation
approximation holds at $T=T_{c}$.

The resummation of the leading order infrared divergences is equivalent with
replacing tree-level bosonic masses in the thermal and Colman-Weinberg parts
of the effective potential with the effective masses Parwani:1991gq

$m^{2}_{i}(h_{I})\rightarrow\overline{m}_{i}^{2}\equiv
m^{2}_{i}(h_{I})+\Pi(T).$ (35)

In case where mass matrix is non-diagonal, the thermal contributions have to
be added prior to mass diagonalization. For instance, the thermally corrected
mass matrix for SM and TS Higgses reads

$\overline{M}_{h,\;H}=\begin{pmatrix}3(\lambda+\kappa+\rho)h_{A}^{2}+\lambda
h_{B}^{2}-f_{0}^{2}(\lambda-\sigma)+\Pi_{h_{A}}&2\lambda h_{A}h_{B}\\\
2\lambda h_{A}h_{B}&3(\lambda+\kappa)h_{B}^{2}+\lambda
h_{A}^{2}-f_{0}^{2}\lambda+\Pi_{h_{B}}\end{pmatrix},$ (36)

while the $Z-\gamma$ mass matrix is

$\overline{M}_{Z,\;\gamma}=\begin{pmatrix}g^{2}h_{A}^{2}+\Pi_{Z}&gg^{\prime}h_{A}^{2}\\\
gg^{\prime}h_{A}^{2}&g^{\prime 2}h_{A}^{2}+\Pi_{\gamma}\end{pmatrix}.$ (37)

Below, we list all thermal corrections $\Pi(T)$ which were used in the
numerical computations. Since these corrections generally depend on the
particle content of the theory, the additional contributions appear for
particular UV extensions considered in this paper. For the ordinary TH the
relevant contributions are

$\displaystyle{\accentset{(\\!\trimbox{0pt1.1ex}{^}\\!)}{\Pi}}_{W}=\frac{11}{6}g{\,}^{2}T^{2},\qquad{\accentset{(\\!\trimbox{0pt1.1ex}{^}\\!)}{\Pi}}_{Z}=\frac{11}{6}g^{2}T^{2},\qquad{\accentset{(\\!\trimbox{0pt1.1ex}{^}\\!)}{\Pi}}_{\gamma}=\frac{11}{6}g^{\prime}{}^{2}T^{2},$
(38)
$\displaystyle\Pi_{h_{A}}=\Big{(}\frac{5}{6}\lambda+\frac{1}{2}\kappa+\frac{1}{2}\rho+\frac{1}{16}g^{\prime
2}+\frac{3}{16}g^{2}+\frac{1}{4}y_{t}^{2}\Big{)}\,T^{2},$
$\displaystyle\Pi_{h_{B}}=\Big{(}\frac{5}{6}\lambda+\frac{1}{2}\kappa+\frac{1}{16}g^{\prime}{}^{2}+\frac{3}{16}g^{2}+\frac{1}{4}\hat{y}_{t}^{2}\Big{)}\,T^{2},$

As in the main text, we use here hat to denote TS particles and couplings. The
hats are suppressed only in the case of twin Higgs, $g$ and $g^{\prime}$ for
which they are redundant. We stress that only the longitudinal degrees of
freedom of the massive vectors receive thermal contributions at the leading
order. The transverse modes receive thermal corrections of order
$\mathcal{O}(g^{4})$ from non-perturbative effects. Those sub-leading non-
perturbative corrections were neglected in our numerical computations.

In the UV agnostic TH model with enhanced $\hat{y}_{l}$, one has to include
leptonic thermal contribution to the twin Higgs mass

$\Pi_{h_{B}}=\Big{(}\frac{5}{6}\lambda+\frac{1}{2}\kappa+\frac{1}{16}g^{\prime}{}^{2}+\frac{3}{16}g^{2}+\frac{1}{4}\hat{y}_{t}^{2}+\frac{\hat{n}_{l}}{48}\hat{y}_{l}^{2}\Big{)}T^{2},$
(39)

where $\hat{n}_{l}$ is the number of degrees of freedom coupled with the
enhanced Yukawa $\hat{y}_{l}$.

As far as light sleptons are considered, one needs to account for their
thermal masses Comelli:1996vm

$\displaystyle\hat{\Pi}_{\text{sl}\;R}=\Big{\\{}\frac{1}{12}g^{\prime
2}\big{[}4+\frac{1}{2}(\hat{n}_{\text{sl, R}}-\hat{n}_{\text{sl,
L}})-\cos(2\beta)\big{]}+\frac{1}{6}\tilde{y}_{l}^{2}(\cos^{2}\beta+\theta_{\tilde{e}_{L}})\Big{\\}}T^{2},$
(40)
$\displaystyle\hat{\Pi}_{\text{sl}\;L}=\Big{\\{}\frac{1}{4}g^{2}+\frac{1}{24}g^{\prime
2}\big{[}2+\frac{1}{2}(\hat{n}_{\text{sl, L}}-\hat{n}_{\text{sl,
R}})+\cos(2\beta)\big{]}+\frac{1}{12}\tilde{y}_{l}^{2}(\cos^{2}\beta+\theta_{\tilde{e}_{R}})\Big{\\}}T^{2}.$

Here, $\hat{n}_{\text{sl, L}}$ ($\hat{n}_{\text{sl, R}}$) is the number of the
left-handed (right-handed) twin sleptonic degrees of freedom coupled to twin
Higgs via enhanced Yukawa coupling. $\theta_{\tilde{e}_{I}}$ is unity if
$\tilde{e}_{I}$ is present in thermal equilibrium and zero if it is decoupled.
In the MSSM-like extension also the thermal masses of TS Higgs and massive
vectors coupled with sleptons are modified Comelli:1996vm

$\displaystyle\Pi_{h_{B}}=\Big{[}\frac{5}{6}\lambda+\frac{1}{2}\kappa+\frac{1}{16}g^{\prime}{}^{2}+\frac{3}{16}g^{2}+\frac{1}{4}\tilde{y}_{t}^{2}+\frac{\hat{n}_{l}}{48}\tilde{y}_{l}^{2}\cos^{2}(\beta)+\frac{\hat{n}_{\text{sl}}}{24}\tilde{y}^{2}_{l}\cos^{2}(\beta)\Big{]}T^{2},$
(41)
$\displaystyle{\accentset{(\\!\trimbox{0pt1.1ex}{^}\\!)}{\Pi}}_{W}={\accentset{(\\!\trimbox{0pt1.1ex}{^}\\!)}{\Pi}}_{Z}=\Big{(}\frac{11}{6}+\frac{1}{12}{\accentset{(\\!\trimbox{0pt1.1ex}{^}\\!)}{n}}_{\text{sl,
L}}\Big{)}g^{2}T^{2},\qquad\quad{\accentset{(\\!\trimbox{0pt1.1ex}{^}\\!)}{\Pi}}_{\gamma}=\Big{(}\frac{11}{6}+\frac{1}{12}{\accentset{(\\!\trimbox{0pt1.1ex}{^}\\!)}{n}}_{\text{sl,
L}}+\frac{1}{6}{\accentset{(\\!\trimbox{0pt1.1ex}{^}\\!)}{n}}_{\text{sl,
R}}\Big{)}g^{\prime}{}^{2}T^{2}.$

where $\hat{n}_{\text{sl}}$ is the total number of twin sleptonic degrees of
freedom coupled to twin Higgs via enhanced Yukawa coupling.

##  Appendix C Estimation of GW signal

###  C.1 Input parameters

In order to obtain the power spectrum of GW, one first needs to compute the
percolation temperature $T_{p}\approx T_{n}$, latent heat of the transition
$\alpha$, (inverse) time duration of the transition in Hubble units
$\beta_{\ast}/H_{\ast}$ and terminal wall velocity $v_{w}$. The first three
parameters can be relatively easy obtained from the effective potential (5)
and three-dimensional Euclidean action $S_{3}$ Lewicki:2021pgr ;
LISACosmologyWorkingGroup:2022jok

$\displaystyle\alpha=\frac{30}{\pi^{2}T_{n}^{4}g_{\ast}}\left(\Delta
V_{\text{eff}}(T_{n})-\frac{1}{4}T_{n}\Delta\left.\frac{\partial
V_{\text{eff}}}{\partial T}\right|_{T=T_{n}}\right),$ (42)
$\displaystyle\frac{\beta_{\ast}}{H_{\ast}}=T\left.\frac{\partial}{\partial
T}\left(\frac{S_{3}}{T}\right)\right|_{T=T_{n}},$ (43)

where $\Delta$ denotes the difference between quantities evaluated in the
false and true vacua while $g_{\ast}$ is the number of relativistic degrees of
freedom in the model at $T=T_{n}$. The terminal bubble wall velocity was
obtained using equilibrium approximation Lewicki:2021pgr

$v_{w}^{2}\approx\frac{\Delta V_{\text{eff}}}{\Delta
V_{\text{eff}}(T_{n})-\frac{1}{4}T_{n}\Delta\left.\frac{\partial
V_{\text{eff}}}{\partial T}\right|_{T=T_{n}}},$ (44)

which turns out to be quite accurate and model independent as long as bubbles
remain subsonic. When the expression (44) exceeds the speed of sound
$c_{s}^{2}=1/3$, one can assume that $v_{w}=1$.

###  C.2 Spectrum

During FOPT there are three main sources of GW: bubble wall collisions,
magnetohydrodynamic turbulence and sound waves produced in bubble collisions.
Their contribution to the overall GW signal highly depends on the phase
transition details.

At all chosen test points, the friction exerted by plasma was sufficient to
stop the acceleration of bubble walls at equilibrium terminal velocity since
the run-away condition Espinosa:2010hh ; Caprini:2015zlo

$\alpha>\alpha_{\infty}\approx\frac{30}{24\pi^{2}}\frac{\sum_{\text{bos}}\Delta
m_{i}^{2}(T=0)+\frac{1}{2}\sum_{\text{ferm}}\Delta
m_{i}^{2}(T=0)}{g_{\ast}T_{n}^{2}}$ (45)

was never fulfilled. Here, $\Delta$ denotes difference between the mass
squared evaluated in the new and the old phase. On the other hand, at all test
points equilibrium terminal velocities obtained with (44) were subsonic. Thus,
we always stay in the regime where bubble wall collisions are a negligible
source, and GW generation is dominated by sound waves Ellis:2018mja . In this
case, the approximate analytic expression for the GW spectrum reads
Caprini:2015zlo

$h^{2}\Omega_{\text{SW}}=2.65\times
10^{-6}\left(\frac{H_{\ast}}{\beta_{\ast}}\right)\left(\frac{\kappa_{\text{SW}}\alpha}{\alpha+1}\right)^{2}\left(\frac{100}{g_{\ast}}\right)^{\frac{1}{3}}v_{w}\,S_{\text{SW}}(f).$
(46)

Here, $\kappa_{\text{SW}}$ is the fraction of the latent energy released in
the transition that was transmitted to sound waves. To find its true value,
one would need to solve the equations of state for the scalar phase. In this
work, we use approximate formulae for $\kappa_{\text{SW}}$ derived in
Espinosa:2010hh  which turn out to be accurate for a large class of models
with polynomial scalar potentials Caprini:2019egz . The profile function for
the sound waves $S_{\text{SW}}(f)$ is given by Caprini:2015zlo

$S_{\text{SW}}(f)=\left(\frac{f}{f_{\text{SW}}}\right)^{3}\bigg{(}\frac{7}{4+3\left(\frac{f}{f_{\text{SW}}}\right)^{2}}\bigg{)}^{\frac{7}{2}},$

where

$f_{\text{SW}}=1.9\times
10^{-5}\text{Hz}\,\frac{1}{v_{w}}\left(\frac{\beta_{\ast}}{H_{\ast}}\right)\left(\frac{T_{\ast}}{100\,\text{GeV}}\right)\left(\frac{g_{\ast}}{100}\right)^{\frac{1}{6}}.$

## References

  * (1) Z. Chacko, H.-S. Goh, and R. Harnik, The Twin Higgs: Natural electroweak breaking from mirror symmetry, Phys. Rev. Lett. 96 (2006) 231802, [hep-ph/0506256].
  * (2) Z. Chacko, Y. Nomura, M. Papucci, and G. Perez, Natural little hierarchy from a partially goldstone twin Higgs, JHEP 01 (2006) 126, [hep-ph/0510273].
  * (3) Z. Chacko, H.-S. Goh, and R. Harnik, A Twin Higgs model from left-right symmetry, JHEP 01 (2006) 108, [hep-ph/0512088].
  * (4) A. Falkowski, S. Pokorski, and M. Schmaltz, Twin SUSY, Phys. Rev. D 74 (2006) 035003, [hep-ph/0604066].
  * (5) S. Chang, L. J. Hall, and N. Weiner, A Supersymmetric twin Higgs, Phys. Rev. D 75 (2007) 035009, [hep-ph/0604076].
  * (6) N. Craig and K. Howe, Doubling down on naturalness with a supersymmetric twin Higgs, JHEP 03 (2014) 140, [arXiv:1312.1341].
  * (7) A. Katz, A. Mariotti, S. Pokorski, D. Redigolo, and R. Ziegler, SUSY Meets Her Twin, JHEP 01 (2017) 142, [arXiv:1611.08615].
  * (8) M. Badziak and K. Harigaya, Supersymmetric D-term Twin Higgs, JHEP 06 (2017) 065, [arXiv:1703.02122].
  * (9) M. Badziak and K. Harigaya, Minimal Non-Abelian Supersymmetric Twin Higgs, JHEP 10 (2017) 109, [arXiv:1707.09071].
  * (10) M. Badziak and K. Harigaya, Asymptotically Free Natural Supersymmetric Twin Higgs Model, Phys. Rev. Lett. 120 (2018), no. 21 211803, [arXiv:1711.11040].
  * (11) N. Craig, S. Knapen, and P. Longhi, Neutral Naturalness from Orbifold Higgs Models, Phys. Rev. Lett. 114 (2015), no. 6 061803, [arXiv:1410.6808].
  * (12) N. Craig, S. Knapen, and P. Longhi, The Orbifold Higgs, JHEP 03 (2015) 106, [arXiv:1411.7393].
  * (13) M. Geller and O. Telem, Holographic Twin Higgs Model, Phys. Rev. Lett. 114 (2015) 191801, [arXiv:1411.2974].
  * (14) P. Batra and Z. Chacko, A Composite Twin Higgs Model, Phys. Rev. D 79 (2009) 095012, [arXiv:0811.0394].
  * (15) R. Barbieri, D. Greco, R. Rattazzi, and A. Wulzer, The Composite Twin Higgs scenario, JHEP 08 (2015) 161, [arXiv:1501.07803].
  * (16) M. Low, A. Tesi, and L.-T. Wang, Twin Higgs mechanism and a composite Higgs boson, Phys. Rev. D 91 (2015) 095012, [arXiv:1501.07890].
  * (17) C. Csaki, M. Geller, O. Telem, and A. Weiler, The Flavor of the Composite Twin Higgs, JHEP 09 (2016) 146, [arXiv:1512.03427].
  * (18) R. Contino, D. Greco, R. Mahbubani, R. Rattazzi, and R. Torre, Precision Tests and Fine Tuning in Twin Higgs Models, Phys. Rev. D 96 (2017), no. 9 095036, [arXiv:1702.00797].
  * (19) R. Barbieri, T. Gregoire, and L. J. Hall, Mirror world at the large hadron collider, hep-ph/0509242.
  * (20) N. Craig, A. Katz, M. Strassler, and R. Sundrum, Naturalness in the Dark at the LHC, JHEP 07 (2015) 105, [arXiv:1501.05310].
  * (21) D. Buttazzo, F. Sala, and A. Tesi, Singlet-like Higgs bosons at present and future colliders, JHEP 11 (2015) 158, [arXiv:1505.05488].
  * (22) A. Ahmed, Heavy Higgs of the Twin Higgs Models, JHEP 02 (2018) 048, [arXiv:1711.03107].
  * (23) I. Garcia Garcia, R. Lasenby, and J. March-Russell, Twin Higgs WIMP Dark Matter, Phys. Rev. D 92 (2015), no. 5 055034, [arXiv:1505.07109].
  * (24) M. Farina, Asymmetric Twin Dark Matter, JCAP 11 (2015) 017, [arXiv:1506.03520].
  * (25) N. Craig and A. Katz, The Fraternal WIMP Miracle, JCAP 10 (2015) 054, [arXiv:1505.07113].
  * (26) R. Barbieri, L. J. Hall, and K. Harigaya, Minimal Mirror Twin Higgs, JHEP 11 (2016) 172, [arXiv:1609.05589].
  * (27) R. Barbieri, L. J. Hall, and K. Harigaya, Effective Theory of Flavor for Minimal Mirror Twin Higgs, JHEP 10 (2017) 015, [arXiv:1706.05548].
  * (28) M. Badziak, G. Grilli Di Cortona, and K. Harigaya, Natural Twin Neutralino Dark Matter, Phys. Rev. Lett. 124 (2020), no. 12 121803, [arXiv:1911.03481].
  * (29) J. Terning, C. B. Verhaaren, and K. Zora, Composite Twin Dark Matter, Phys. Rev. D 99 (2019), no. 9 095020, [arXiv:1902.08211].
  * (30) Z. Chacko, D. Curtin, M. Geller, and Y. Tsai, Direct detection of mirror matter in Twin Higgs models, JHEP 11 (2021) 198, [arXiv:2104.02074].
  * (31) M. Badziak, G. Grilli di Cortona, K. Harigaya, and M. Łukawski, Charged dark matter in supersymmetric twin Higgs models, JHEP 10 (2022) 057, [arXiv:2202.10488].
  * (32) V. Prilepina and Y. Tsai, Reconciling Large And Small-Scale Structure In Twin Higgs Models, JHEP 09 (2017) 033, [arXiv:1611.05879].
  * (33) Z. Chacko, D. Curtin, M. Geller, and Y. Tsai, Cosmological Signatures of a Mirror Twin Higgs, JHEP 09 (2018) 163, [arXiv:1803.03263].
  * (34) S. Bansal, J. H. Kim, C. Kolda, M. Low, and Y. Tsai, Mirror twin Higgs cosmology: constraints and a possible resolution to the H0 and S8 tensions, JHEP 05 (2022) 050, [arXiv:2110.04317].
  * (35) Z. Chacko, N. Craig, P. J. Fox, and R. Harnik, Cosmology in Mirror Twin Higgs and Neutrino Masses, JHEP 07 (2017) 023, [arXiv:1611.07975].
  * (36) N. Craig, S. Koren, and T. Trott, Cosmological Signals of a Mirror Twin Higgs, JHEP 05 (2017) 038, [arXiv:1611.07977].
  * (37) K. Harigaya, R. Mcgehee, H. Murayama, and K. Schutz, A predictive mirror twin Higgs with small Z2 breaking, JHEP 05 (2020) 155, [arXiv:1905.08798].
  * (38) M. Kamionkowski, A. Kosowsky, and M. S. Turner, Gravitational radiation from first order phase transitions, Phys. Rev. D 49 (1994) 2837–2851, [astro-ph/9310044].
  * (39) J. Ellis, M. Lewicki, and J. M. No, On the Maximal Strength of a First-Order Electroweak Phase Transition and its Gravitational Wave Signal, JCAP 04 (2019) 003, [arXiv:1809.08242].
  * (40) C. Caprini et al., Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions, JCAP 04 (2016) 001, [arXiv:1512.06239].
  * (41) C. Caprini et al., Detecting gravitational waves from cosmological phase transitions with LISA: an update, JCAP 03 (2020) 024, [arXiv:1910.13125].
  * (42) V. A. Kuzmin, V. A. Rubakov, and M. E. Shaposhnikov, On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe, Phys. Lett. B 155 (1985) 36.
  * (43) M. E. Shaposhnikov, Possible Appearance of the Baryon Asymmetry of the Universe in an Electroweak Theory, JETP Lett. 44 (1986) 465–468.
  * (44) J. M. Cline, Baryogenesis, in Les Houches Summer School - Session 86: Particle Physics and Cosmology: The Fabric of Spacetime, 9, 2006\.  hep-ph/0609145.
  * (45) D. E. Morrissey and M. J. Ramsey-Musolf, Electroweak baryogenesis, New J. Phys. 14 (2012) 125003, [arXiv:1206.2942].
  * (46) P. B. Arnold and O. Espinosa, The Effective potential and first order phase transitions: Beyond leading-order, Phys. Rev. D 47 (1993) 3546, [hep-ph/9212235]. [Erratum: Phys.Rev.D 50, 6662 (1994)].
  * (47) K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov, The Electroweak phase transition: A Nonperturbative analysis, Nucl. Phys. B 466 (1996) 189–258, [hep-lat/9510020].
  * (48) K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov, Is there a hot electroweak phase transition at $m_{H}\gtrsim m_{W}$?, Phys. Rev. Lett. 77 (1996) 2887–2890, [hep-ph/9605288].
  * (49) K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov, A Nonperturbative analysis of the finite T phase transition in SU(2) x U(1) electroweak theory, Nucl. Phys. B 493 (1997) 413–438, [hep-lat/9612006].
  * (50) V. Barger, P. Langacker, M. McCaskey, M. J. Ramsey-Musolf, and G. Shaughnessy, LHC Phenomenology of an Extended Standard Model with a Real Scalar Singlet, Phys. Rev. D 77 (2008) 035005, [arXiv:0706.4311].
  * (51) T. Brauner, T. V. I. Tenkanen, A. Tranberg, A. Vuorinen, and D. J. Weir, Dimensional reduction of the Standard Model coupled to a new singlet scalar field, JHEP 03 (2017) 007, [arXiv:1609.06230].
  * (52) D. Curtin, P. Meade, and C.-T. Yu, Testing Electroweak Baryogenesis with Future Colliders, JHEP 11 (2014) 127, [arXiv:1409.0005].
  * (53) K. Fujikura, K. Kamada, Y. Nakai, and M. Yamaguchi, Phase Transitions in Twin Higgs Models, JHEP 12 (2018) 018, [arXiv:1810.00574].
  * (54) H. Beauchesne, K. Earl, and T. Grégoire, The spontaneous $\mathbb{Z}_{2}$ breaking Twin Higgs, JHEP 01 (2016) 130, [arXiv:1510.06069].
  * (55) S. R. Coleman and E. J. Weinberg, Radiative Corrections as the Origin of Spontaneous Symmetry Breaking, Phys. Rev. D 7 (1973) 1888–1910.
  * (56) M. Quiros, Finite temperature field theory and phase transitions, in ICTP Summer School in High-Energy Physics and Cosmology, pp. 187–259, 1, 1999.  hep-ph/9901312.
  * (57) R. R. Parwani, Resummation in a hot scalar field theory, Phys. Rev. D 45 (1992) 4695, [hep-ph/9204216]. [Erratum: Phys.Rev.D 48, 5965 (1993)].
  * (58) C. Kilic and S. Swaminathan, Can A Pseudo-Nambu-Goldstone Higgs Lead To Symmetry Non-Restoration?, JHEP 01 (2016) 002, [arXiv:1508.05121].
  * (59) P. B. Arnold, The Electroweak phase transition: Part 1. Review of perturbative methods, in 8th International Seminar on High-energy Physics, 10, 1994.  hep-ph/9410294.
  * (60) A. D. Linde, Infrared Problem in Thermodynamics of the Yang-Mills Gas, Phys. Lett. B 96 (1980) 289–292.
  * (61) C. L. Wainwright, CosmoTransitions: Computing Cosmological Phase Transition Temperatures and Bubble Profiles with Multiple Fields, Comput. Phys. Commun. 183 (2012) 2006–2013, [arXiv:1109.4189].
  * (62) A. Beniwal, M. Lewicki, J. D. Wells, M. White, and A. G. Williams, Gravitational wave, collider and dark matter signals from a scalar singlet electroweak baryogenesis, JHEP 08 (2017) 108, [arXiv:1702.06124].
  * (63) P. D. Group, R. Workman, V. Burkert, V. Crede, E. Klempt, U. Thoma, L. Tiator, K. Agashe, G. Aielli, B. Allanach, et al., Review of particle physics, Progress of Theoretical and Experimental Physics 2022 (2022), no. 8 083C01.
  * (64) O. Matsedonskyi and G. Servant, High-Temperature Electroweak Symmetry Non-Restoration from New Fermions and Implications for Baryogenesis, JHEP 09 (2020) 012, [arXiv:2002.05174].
  * (65) O. Matsedonskyi, High-Temperature Electroweak Symmetry Breaking by SM Twins, JHEP 04 (2021) 036, [arXiv:2008.13725].
  * (66) S. P. Martin, A Supersymmetry primer, Adv. Ser. Direct. High Energy Phys. 18 (1998) 1–98, [hep-ph/9709356].
  * (67) ALEPH Collaboration, A. Heister et al., Search for scalar leptons in e+ e- collisions at center-of-mass energies up to 209-GeV, Phys. Lett. B 526 (2002) 206–220, [hep-ex/0112011].
  * (68) DELPHI Collaboration, J. Abdallah et al., Searches for supersymmetric particles in e+ e- collisions up to 208-GeV and interpretation of the results within the MSSM, Eur. Phys. J. C 31 (2003) 421–479, [hep-ex/0311019].
  * (69) OPAL Collaboration, G. Abbiendi et al., Search for anomalous production of dilepton events with missing transverse momentum in e+ e- collisions at s**(1/2) = 183-Gev to 209-GeV, Eur. Phys. J. C 32 (2004) 453–473, [hep-ex/0309014].
  * (70) L3 Collaboration, P. Achard et al., Search for scalar leptons and scalar quarks at LEP, Phys. Lett. B 580 (2004) 37–49, [hep-ex/0310007].
  * (71) CMS Collaboration, Search for direct pair production of supersymmetric partners of $\tau$ leptons in the final state with two hadronically decaying $\tau$ leptons and missing transverse momentum in proton-proton collisions at $\sqrt{s}$ = 13 TeV, arXiv:2207.02254.
  * (72) ATLAS Collaboration, G. Aad et al., Search for electroweak production of charginos and sleptons decaying into final states with two leptons and missing transverse momentum in $\sqrt{s}=13$ TeV $pp$ collisions using the ATLAS detector, Eur. Phys. J. C 80 (2020), no. 2 123, [arXiv:1908.08215].
  * (73) ATLAS Collaboration, G. Aad et al., Searches for electroweak production of supersymmetric particles with compressed mass spectra in $\sqrt{s}=$ 13 TeV $pp$ collisions with the ATLAS detector, Phys. Rev. D 101 (2020), no. 5 052005, [arXiv:1911.12606].
  * (74) ATLAS Collaboration, Search for direct pair production of sleptons and charginos decaying to two leptons and neutralinos with mass splittings near the $W$ boson mass in ${\sqrt{s}=13\,}$TeV $pp$ collisions with the ATLAS detector, .
  * (75) X. Wang, F. P. Huang, and X. Zhang, Phase transition dynamics and gravitational wave spectra of strong first-order phase transition in supercooled universe, JCAP 05 (2020) 045, [arXiv:2003.08892].
  * (76) M. Merchand, private communication.
  * (77) A. Bertoldi et al., AEDGE: Atomic experiment for dark matter and gravity exploration in space, Exper. Astron. 51 (2021), no. 3 1417–1426, [arXiv:1908.00802].
  * (78) M. Punturo et al., The Einstein Telescope: A third-generation gravitational wave observatory, Class. Quant. Grav. 27 (2010) 194002\.
  * (79) S. Hild et al., Sensitivity Studies for Third-Generation Gravitational Wave Observatories, Class. Quant. Grav. 28 (2011) 094013, [arXiv:1012.0908].
  * (80) E. Thrane and J. D. Romano, Sensitivity curves for searches for gravitational-wave backgrounds, Phys. Rev. D 88 (2013), no. 12 124032, [arXiv:1310.5300].
  * (81) L. S. Collaboration et al., Advanced ligo, Class. Quantum Gravity 32 (2015), no. 7 074001, [arXiv:1411.4547].
  * (82) N. Bartolo et al., Science with the space-based interferometer LISA. IV: Probing inflation with gravitational waves, JCAP 12 (2016) 026, [arXiv:1610.06481].
  * (83) C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni, A. Ricciardone, and G. Tasinato, Reconstructing the spectral shape of a stochastic gravitational wave background with LISA, JCAP 11 (2019) 017, [arXiv:1906.09244].
  * (84) L. Badurina et al., AION: An Atom Interferometer Observatory and Network, JCAP 05 (2020) 011, [arXiv:1911.11755].
  * (85) M. Lewicki, private communication.
  * (86) D. Curtin, P. Meade, and H. Ramani, Thermal Resummation and Phase Transitions, Eur. Phys. J. C 78 (2018), no. 9 787, [arXiv:1612.00466].
  * (87) D. Comelli and J. R. Espinosa, Bosonic thermal masses in supersymmetry, Phys. Rev. D 55 (1997) 6253–6263, [hep-ph/9606438].
  * (88) 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].
  * (89) LISA Cosmology Working Group Collaboration, P. Auclair et al., Cosmology with the Laser Interferometer Space Antenna, arXiv:2204.05434.
  * (90) J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, Energy Budget of Cosmological First-order Phase Transitions, JCAP 06 (2010) 028, [arXiv:1004.4187].