# Probing intermediate scale Froggatt-Nielsen models at future gravitational
wave observatories

Dhruv Ringe<EMAIL_ADDRESS>Discipline of Physics, IIT Indore,
Khandwa Road, Simrol, Indore, Madhya Pradesh, India -453552

###### Abstract

The flavor symmetry-breaking scale in the Froggatt-Nielsen (FN) mechanism is
very weakly constrained by present experiments, and can lie anywhere from a
few TeV to the Planck scale. We construct two minimal, non-supersymmetric,
ultraviolet (UV) complete models that generate the FN mechanism, with a global
$U(1)_{\rm{FN}}$ flavor symmetry and a single flavon field. Using the one-loop
finite temperature effective potential, we explore the possibility of a strong
first order phase transition (SFOPT) induced by the flavon. We show that if
the flavor symmetry-breaking occurs at intermediate scales $\sim
10^{4}-10^{7}$ GeV, then in certain regions of the parameter space, the
associated stochastic gravitational wave (GW) background is strong enough to
be detected by second generation GW observatories such as the Big Bang
Observer (BBO), the Deci-hertz interferometer Gravitational Observatory
(DECIGO), the Cosmic Explorer (CE) and the Einstein Telescope (ET). We
identify viable regions of the parameter space for the best detection
prospects. While both models of flavor can produce a detectable GW background,
the GW signature cannot be used to discriminate between them.

##  1 Introduction

The observation of the first GW signal from a binary neutron star merger
GW2015 [1] by the aLIGO [2] collaboration has opened up a new window to
explore the early history of the universe. Our capability of observing feeble
GW signals is expected to be significantly enhanced by several planned GW
observatories such as LISA [3], BBO [4], DECIGO [5], CE [6], and ET [7]. Along
with pulsar timing arrays (PTA’s) [8, 9, 10, 11, 12, 13], these observatories
will detect GW signals with frequencies ranging from $\sim 10^{-9}$ Hz to
$\sim 10^{4}$ Hz. With binary neurton star and binary black hole merger events
routinely being detected by LIGO and VIRGO collaborations [14, 15, 16], the
prospect of detecting a stochastic cosmological background of GW’s by these
new obervatories is indeed promising. Several possible sources of a stochastic
GW background have been proposed, including the density fluctuations from
inflation [17, 18, 19, 20, 21], topological defects in the vacuum [22, 23, 24,
25, 26, 27], and cosmological phase transitions [28, 29].

The final piece of the minimal Standard Model (SM) of particle physics was
found with the discovery of the Higgs boson at the Large Hadron Collider (LHC)
in 2012 [30, 31]. Despite its tremendous success, the SM has limitations which
hint at possible new physics. The SM does not explain the hierarchy of fermion
masses, or equivalently the hierarchy of yukawa couplings, which range from
$Y\sim 1$ for the top-quark, to $Y\sim 10^{-6}$ for the electron. Accounting
for neutrino masses further enhances this hierarchy, making the problem much
worse. Equally puzzling is the hierarchy of quark mixing angles, and the
anarchy of leptonic mixing angles. This problem of yukawa hierarchies has come
to be known as the SM flavor puzzle [32, 33, 34].

In this paper, we revisit the elegant solution to the SM flavor puzzle,
originially proposed by Froggatt and Nielsen [35], focussing on the quark
sector. The observed hiearchy of quark masses and mixings is generated by non-
renormalizable operators involving one or more complex scalar fields
(flavons). The flavons spontaneously break the flavor symmetry imposed on the
SM fields by acquiring a non-zero vacuum expectation value (VEV). A crucial
parameter in FN models is the ratio $\epsilon=\frac{\langle
S\rangle}{\Lambda}$, where $\langle S\rangle$ is the VEV of the flavon, and
$\Lambda$ is the cutoff scale of the effective theory, associated with the
mass of the messenger vector-like fermions (VLF’s) in the UV completion.

The UV completion of FN is a tedious task due to the involvement of several
VLF’s, and is further complicated for realizations of FN with large
$U(1)_{\rm{FN}}$ charge assignments. Some explicit UV constructions are found
for example in [36, 37, 38, 39]. While the method to determine the minimum
number vectorlike fermions was first given in [36], a simple procedure was
given in [38] to find the complete set of charge assignments. Using the
procedure of [38], we construct minimal, non-supersymmetric, UV complete
models for two simple realizations of the FN mechanism with one flavon. It
should be noted however that adding several heavy fermions threatens the
vacuum stability of the flavon potential, which demands the inclusion of
additional bosons to restore the vacuum stability.

Given an FN model, the ratio $\epsilon$ is fixed to some extent by the
observed pattern of quark masses and mixing angles, whereas the scale
$\Lambda$ itself is undetermined. Phenomenology of neutral meson mixing puts a
lower bound of a few TeV on $\Lambda$ [40]. On the other hand, searches for
vector-like fermions based on collider results [41, 42] and electroweak vacuum
stability [43, 44, 45] put a bound of $\gtrsim 1$ TeV. Planned future
colliders such as the International Linear Collider (ILC) [46, 47] and the
Future Circular Collider (FCC) [48, 49, 50] may push the lower bound on
$\Lambda$ to $\sim 50$ TeV. Beyond that, the cutoff is practically
unconstrained and could lie anywhere between $\sim 10$ TeV and the Planck
scale $(\sim 10^{19}$ GeV). However, if the flavon undergoes a SFOPT, then the
corresponding stochastic background may be detected by upcoming GW detectors,
for symmetry-breaking scales as high as $10^{4}-10^{7}$ GeV. The corresponding
GW signature has a peak frequency of $10^{-1}-10^{2}$ Hz, making it an ideal
energy scale to be probed at observatiories like BBO, DECIGO, ET and CE. While
several authors have discussed the GW background arising at such intermediate
scales in different contexts [51, 52, 53, 54, 55, 56, 57, 58, 59], as far as
we know, this is the first attempt to probe FN models incorporating UV
completions, with GW’s.

This paper is organized as follows: in section 2, we review the FN mechanism
and discuss two realizations involving a single flavon. In section 3, we
construct the UV theory for the two FN models by giving the particle content
involving heavy vector-like quarks. In section 4, we construct the one-loop
finite temperature effective potential needed to analyse the flavon phase
transition (PT). The issue of vacuum stability of FN UV completions is
discussed in section 5. In section 6 we define the parameters characterizing a
first order phase transition (FOPT), and give the parameter scans showing
regions of SFOPT in the two models considered. In section 7, we discuss the
prospects of dectecting the GW signature at the future GW observatories.
Finally, we give our concluding remarks in section 8.

##  2 Froggatt-Nielsen mechanism

In its simplest form, the FN mechanism is achieved by imposing a global
horizontal $U(1)_{\rm{FN}}$ symmetry on the SM fields. This symmetry prohibits
the inclusion of the usual yukawa terms in the SM Lagrangian. Instead, we can
write $U(1)_{\rm{FN}}$ preserving effective operators of the form:

$-\mathcal{L}_{\rm{FN}}=y_{ij}\left(\frac{S}{\Lambda}\right)^{n_{ij}}\overline{f}_{Li}Hf_{Rj}+\rm{h.c.}$
(1)

where $y_{ij}$ are $\mathcal{O}(1)$ couplings, $f_{L}$, $f_{R}$ respectively
represent left handed fermion doublets and right handed fermion singlets of
the SM, and $i,j$ denote the family index. $\Lambda$ is the cutoff scale,
related to the mass scale of particles of the UV theory. $S$ is a complex
scalar field called the flavon, and $H$ is the SM Higgs doublet. The FN
mechanism works both in the lepton sector as well as the quark sector. Here we
focus exlusively on the quark sector, where,

$-\mathcal{L}_{\rm{FN}}=y^{u}_{ij}\left(\frac{S}{\Lambda}\right)^{n^{u}_{ij}}\overline{Q}^{i}_{L}\widetilde{H}~{}u^{j}_{R}+y^{d}_{ij}\left(\frac{S}{\Lambda}\right)^{n^{d}_{ij}}\overline{Q}^{i}_{L}H~{}d^{~{}j}_{R}+\rm{h.c.}$
(2)

where $\widetilde{H}=i\sigma^{2}H$. Without loss of generality, we can choose
the FN charges of the scalars as: $Q_{\rm{FN}}(S)=-1$ and $Q_{\rm{FN}}(H)=0$.
The power of the $\frac{S}{\Lambda}$ term is:

$n^{u/d}_{ij}=Q_{\rm{FN}}(\overline{Q}_{L}^{i})+Q_{\rm{FN}}(u_{R}^{~{}j}/d_{R}^{~{}j}).$

When $U(1)_{\rm{FN}}$ is spontaneously broken by the non-zero VEV of the
flavon, the FN operators give rise to the familiar SM yukawa terms:

$-\mathcal{L}_{\rm{Y}}=Y^{u}_{ij}\overline{Q}^{i}_{L}\widetilde{H}~{}u^{j}_{R}+Y^{d}_{ij}\overline{Q}^{i}_{L}H~{}d^{~{}j}_{R}+\rm{h.c.},$
(3)

where ,

$Y^{u/d}_{ij}=y^{u/d}_{ij}\epsilon^{n^{u/d}_{ij}},$ (4)

and $\epsilon$ is defined as the ratio,

$\epsilon\equiv\frac{\langle S\rangle}{\Lambda}=\frac{v_{s}}{\sqrt{2}\Lambda}$
(5)

For $\epsilon<1$, the desired hierarchy of quark masses is generated with
appropriate charge assigments. Corresponding quark mass matrices are:

${\bf{\mathcal{M}^{u/d}}}=\frac{v}{\sqrt{2}}{\bf{Y^{u/d}}},$

where $\frac{v}{\sqrt{2}}$ is the electroweak VEV with $v=246.2$ GeV.

It is useful to parameterize the hierarchy in terms of the Cabbibo angle:
$\lambda=\sin\theta_{cb}\sim 0.23$. The effective couplings $y^{u/d}_{ij}$ are
all required to be $\mathcal{O}(1)$, but their exact value is unimportant,
allowing some freedom in assigning charges to the quarks.

Below we consider two simple examples which implement the FN mechanism with a
single flavon.

  * •

Model 1: $\epsilon=\lambda=0.23$, [60]

$Q_{\rm{FN}}(\overline{Q}_{L})=(3,2,0),~{}Q_{\rm{FN}}(u_{R})=(5,2,0),~{}Q_{\rm{FN}}(d_{R})=(4,3,3).$
(6)

The yukawa matrices are:

$\displaystyle{\bf{Y^{u}}}\sim\begin{pmatrix}\epsilon^{8}&\epsilon^{5}&\epsilon^{3}\\\
\epsilon^{7}&\epsilon^{4}&\epsilon^{2}\\\
\epsilon^{5}&\epsilon^{2}&1\end{pmatrix},~{}~{}~{}{\bf{Y^{d}}}\sim\begin{pmatrix}\epsilon^{7}&\epsilon^{6}&\epsilon^{6}\\\
\epsilon^{6}&\epsilon^{5}&\epsilon^{5}\\\
\epsilon^{4}&\epsilon^{3}&\epsilon^{3}\end{pmatrix},$ (7)

where each entry of the matrices is multiplied by an $\mathcal{O}(1)$ factor
$y^{u/d}_{ij}$. It is useful to note the $\epsilon$-dependence of the
determinants:

$\det{\bf{Y^{u}}}\propto\epsilon^{12},~{}~{}~{}\det{\bf{Y^{d}}}\propto\epsilon^{15}.$
(8)

  * •

Model 2: This is a non-supersymmetric variation of the charge assignment
considered in [61]. $\epsilon=\lambda^{2}=0.05$,

$Q_{\rm{FN}}(\overline{Q}_{L})=Q_{\rm{FN}}(u_{R})=(2,1,0),~{}Q_{\rm{FN}}(d_{R})=(1,1,1).$
(9)

The yukawa matrices are 111Note that the matrix $\bf{Y^{d}}$ differs from the
matrix in [61] by a factor of $\epsilon$. This is because of the presence of
two Higgs doublets in supersymmetric theories, which can take care of the
additional hierarchy.:

$\displaystyle{\bf{Y^{u}}}\sim\begin{pmatrix}\epsilon^{4}&\epsilon^{3}&\epsilon^{2}\\\
\epsilon^{3}&\epsilon^{2}&\epsilon\\\
\epsilon^{2}&\epsilon&1\end{pmatrix},~{}~{}~{}{\bf{Y^{d}}}\sim\begin{pmatrix}\epsilon^{3}&\epsilon^{3}&\epsilon^{3}\\\
\epsilon^{2}&\epsilon^{2}&\epsilon^{2}\\\
\epsilon&\epsilon&\epsilon\end{pmatrix}.$ (10)

The determinants are:

$\det{\bf{Y^{u}}}\propto\epsilon^{6},~{}~{}~{}\det{\bf{Y^{d}}}\propto\epsilon^{6}.$
(11)

The FN mechanism can also be realized with more complicated abelian/non-
abelian groups instead of $U(1)_{\rm{FN}}$, that may or may not be gauged. As
a result, there can be models with more than one flavon field. See Appendix A
for an example of a two flavon model. As we’ll see in the next section, lower
powers of $\epsilon$ imply the need for fewer heavy quarks, and hence a
simpler UV theory.

##  3 UV completion of FN models

Figure 1:  A chain diagram of the UV theory, contributing to an FN operator
corresponding to $Y^{u/d}_{ij}$. The number of insertions of the $S$ field is
equal to the power of $\epsilon$ in $Y^{u/d}_{ij}$. Here $F$’s represent the
heavy VLQ’s sitting at the scale $\Lambda$. Yukawa couplings at the vertices
are all $\mathcal{O}(1)$. Arrows indicate the flow of $U(1)_{\rm{FN}}$ charge,
which is conserved at each vertex.

In the high energy theory, the effective operators listed in (2). are
generated from chain diagrams of the kind shown in figure 1. These diagrams
are mediated by a set of heavy vector-like quarks (VLQ’s) with mass
$M\sim\Lambda$. We assume all VLQ’s to have the same bare mass $M$. The heavy
VLQ’s can either be $SU(2)_{L}$ singlets, or doublets. For simplicity we take
all the VLQ’s to be $SU(2)_{L}$ singlets, as it does not impact our ultimate
goal of analysing the flavon phase transition. We represent the chain diagram
of figure 1 using the short hand notation,

$u^{j}/d^{j}-F^{a}-F^{b}-\cdots-F^{m}-F^{n}-\overline{Q}^{i}$

The scalar field at each vertex is suppressed, and only the SM fermions and
the mediating VLQ’s are shown. For $SU(2)_{L}$ singlet VLQ’s, it is assumed
that the scalar attached to $\overline{Q}$ is $H/\widetilde{H}$, while an $S$
field is attached to each of the other vertices.

The singlet VLQ’s come in ‘up-type’ or ‘down-type’ representations of the SM
gauge group, $\mathcal{G}_{\rm{SM}}=SU(3)_{c}\times SU(2)_{L}\times U_{1}(Y),$
depending upon the type of light quark they interact with. The notation used
here is:

$(T_{L},T_{R})\sim({\bf{3},\bf{1}},4/3)~{}~{}({\rm{up\mathchar 45\relax
type}}),$ $(B_{L},B_{R})\sim({\bf{3},\bf{1}},-2/3)~{}~{}({\rm{down\mathchar
45\relax type}}).$

The fermionic mass Lagrangian can be written as,

$-\mathcal{L}_{\rm{UV}}\supset\overline{F^{u}}_{L}\mathcal{M}^{u}F^{u}_{R}+\overline{F^{d}}_{L}\mathcal{M}^{d}F^{d}_{R}+\rm{h.c.},$
(12)

where $F^{u}_{L}=\big{(}Q^{1},Q^{2},Q^{3},T^{1},T^{2}\cdots\big{)}_{L}^{T}$,
$F^{u}_{R}=\big{(}u,c,t,T^{1},T^{2}\cdots\big{)}_{R}^{T}$,
$F^{d}_{L}=\big{(}Q^{1},Q^{2},Q^{3},B^{1},B^{2}\cdots\big{)}_{L}^{T}$,
$F^{d}_{L}=\big{(}d,s,b,B^{1},B^{2}\cdots\big{)}_{R}^{T}$.

The mass matrices $\mathcal{M}^{u/d}$ depend on the specific model. Physical
fermion masses are given by the singular values of $\mathcal{M}^{u/d}$, i.e.
eigenvalues of $\big{(}\mathcal{M}^{\dagger}\mathcal{M}\big{)}^{1/2}$. In [36,
37], it was shown that for a yukawa matrix Y of the SM fermions with
$\det\textbf{Y}\propto\epsilon^{n}$, at least $n$ VLQ’s are required in the UV
theory. We now follow the procedure elucidated in [38] to assign charges to
the VLQ’s, and build the particle content for our UV models.

###  3.1 Model 1

According to (8), atleast $12$ u-type VLQ’s and $15$ d-type VLQ’s are needed.
For the up sector, first we select three elements from the $\textbf{Y}^{u}$
matrix in (7), which contribute to the determinant 222Every term in the
determinant comes with the same power of $\epsilon$. This is a general feature
of choosing a $U(1)$ symmetry.. For each element, we form a chain diagram with
the number of VLQ’s equal to the corresponding power of $\epsilon$. Here we
choose, $Y^{u}_{11}\propto\epsilon^{8}$, $Y^{u}_{22}\propto\epsilon^{4}$, and
$Y^{u}_{33}\propto\epsilon^{0}$. Next, we use a separate set of VLQ’s to
construct each chain. The $U(1)_{\rm{FN}}$ charge assignment is determined by
conserving the charge at at every vertex.

$\displaystyle
u_{R}^{1}-T^{8}-T^{7}-T^{6}-T^{5}-T^{4}-T^{3}-T^{2}-T^{1}-\overline{Q}_{L}^{1}$
$\displaystyle u_{R}^{2}-T^{\prime 4}-T^{\prime 3}-T^{\prime 2}-T^{\prime
1}-\overline{Q}_{L}^{2}$

Similarly, for the down sector, we have: $Y^{d}_{11}\propto\epsilon^{7}$,
$Y^{d}_{22}\propto\epsilon^{5}$, and $Y^{d}_{33}\propto\epsilon^{3}$. The
chains are,

$\displaystyle
d_{R}^{1}-B^{7}-B^{6}-B^{5}-B^{4}-B^{3}-B^{2}-B^{1}-\overline{Q}_{L}^{1}$
$\displaystyle d_{R}^{2}-B^{\prime 5}-B^{\prime 4}-B^{\prime 3}-B^{\prime
2}-B^{\prime 1}-\overline{Q}_{L}^{2}$ $\displaystyle d_{R}^{3}-B^{\prime\prime
3}-B^{\prime\prime 2}-B^{\prime\prime 1}-\overline{Q}_{L}^{3}$

$Q_{\rm{FN}}$ | -3 | -2 | -1 | 0 | 1 | 2 | 3 | 4
---|---|---|---|---|---|---|---|---
up-type | $T^{1}$ | $T^{2}$ | $T^{3}$ | $T^{4}$ | $T^{5}$ | $T^{6}$ | $T^{7}$ | $T^{8}$
|  | $T^{\prime 1}$ | $T^{\prime 2}$ | $T^{\prime 3}$ | $T^{\prime 4}$ |  |  |
down-type | $B^{1}$ | $B^{2}$ | $B^{3}$ | $B^{4}$ | $B^{5}$ | $B^{6}$ | $B^{7}$ |
|  | $B^{\prime 1}$ | $B^{\prime 2}$ | $B^{\prime 3}$ | $B^{\prime 4}$ | $B^{\prime 5}$ |  |
|  |  |  | $B^{\prime\prime 1}$ | $B^{\prime\prime 2}$ | $B^{\prime\prime 3}$ |  |
Table 1: $U(1)_{\rm{FN}}$ charges for VLQ’s in model 1. Total $12+15=27$ heavy
quarks are needed in this model.

Chains for the rest of the terms can now be written using the VLQ’s found
above. The $U(1)_{\rm{FN}}$ charge assignments for the VLQ’s are displayed in
table 1. The resulting mass matrices are,

$\mathcal{M}^{u}=\setcounter{MaxMatrixCols}{15}\begin{pmatrix}0&0&0&\langle\widetilde{H}\rangle&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&\langle\widetilde{H}\rangle&\langle\widetilde{H}\rangle&0&0&0&0&0&0&0&0&0\\\
0&0&\langle\widetilde{H}\rangle&0&0&0&0&0&\langle\widetilde{H}\rangle&\langle\widetilde{H}\rangle&0&0&0&0&0\\\
0&0&0&M&\langle S\rangle&\langle S\rangle&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&M&0&\langle S\rangle&\langle S\rangle&0&0&0&0&0&0&0\\\
0&0&0&0&0&M&\langle S\rangle&\langle S\rangle&0&0&0&0&0&0&0\\\ 0&0&\langle
S\rangle&0&0&0&M&0&\langle S\rangle&\langle S\rangle&0&0&0&0&0\\\ 0&0&\langle
S\rangle&0&0&0&0&M&\langle S\rangle&\langle S\rangle&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&M&0&\langle S\rangle&\langle S\rangle&0&0&0\\\
0&0&0&0&0&0&0&0&0&M&\langle S\rangle&\langle S\rangle&0&0&0\\\ 0&\langle
S\rangle&0&0&0&0&0&0&0&0&M&0&\langle S\rangle&0&0\\\ 0&\langle
S\rangle&0&0&0&0&0&0&0&0&0&M&\langle S\rangle&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&M&\langle S\rangle&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&M&\langle S\rangle\\\ \langle
S\rangle&0&0&0&0&0&0&0&0&0&0&0&0&0&M\\\ \end{pmatrix},$ (13a)
$\mathcal{M}^{d}=\setcounter{MaxMatrixCols}{18}\begin{pmatrix}0&0&0&\langle
H\rangle&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&\langle H\rangle&\langle
H\rangle&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&\langle H\rangle&\langle
H\rangle&\langle H\rangle&0&0&0&0&0&0&0\\\ 0&0&0&M&\langle S\rangle&\langle
S\rangle&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&M&0&\langle S\rangle&\langle
S\rangle&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&M&\langle S\rangle&\langle
S\rangle&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&M&0&\langle S\rangle&\langle
S\rangle&\langle S\rangle&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&M&\langle
S\rangle&\langle S\rangle&\langle S\rangle&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&M&0&0&\langle S\rangle&\langle S\rangle&\langle
S\rangle&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&M&0&\langle S\rangle&\langle
S\rangle&\langle S\rangle&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&M&\langle
S\rangle&\langle S\rangle&\langle S\rangle&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&M&0&0&\langle S\rangle&\langle S\rangle&\langle
S\rangle&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&M&0&\langle S\rangle&\langle
S\rangle&\langle S\rangle&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&M&\langle
S\rangle&\langle S\rangle&\langle S\rangle&0\\\ 0&\langle S\rangle&\langle
S\rangle&0&0&0&0&0&0&0&0&0&0&0&M&0&0&\langle S\rangle\\\ 0&\langle
S\rangle&\langle S\rangle&0&0&0&0&0&0&0&0&0&0&0&0&M&0&\langle S\rangle\\\
0&\langle S\rangle&\langle S\rangle&0&0&0&0&0&0&0&0&0&0&0&0&0&M&\langle
S\rangle\\\ \langle S\rangle&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&M\end{pmatrix}.$
(13b)

The first three rows and columns of each matrix correspond to the respective
light quarks in UV theory. The VLQ’s are arranged in order of increasing
$Q_{\rm{FN}}$.

###  3.2 Model 2

According to (11), atleast $6$ u-type VLQ’s and $6$ d-type VLQ’s are needed.
Repeating the procedure described above, for the up sector we have:
$Y^{u}_{11}\propto\epsilon^{4}$, $Y^{u}_{22}\propto\epsilon^{2}$, and
$Y^{u}_{33}\propto\epsilon^{0}$. Therefore we have the chains,

$\displaystyle u_{R}^{1}-T^{4}-T^{3}-T^{2}-T^{1}-\overline{Q}_{L}^{1}$
$\displaystyle u_{R}^{2}-T^{\prime 2}-T^{\prime 1}-\overline{Q}_{L}^{2}$

For the down sector, we have: $Y^{d}_{11}\propto\epsilon^{3}$,
$Y^{d}_{22}\propto\epsilon^{2}$, and $Y^{d}_{33}\propto\epsilon^{1}$. This
gives,

$\displaystyle d_{R}^{1}-B^{3}-B^{2}-B^{1}-\overline{Q}_{L}^{1}$
$\displaystyle d_{R}^{2}-B^{\prime 2}-B^{\prime 1}-\overline{Q}_{L}^{2}$
$\displaystyle d_{R}^{3}-B^{\prime\prime 1}-\overline{Q}_{L}^{3}$

$Q_{\rm{FN}}$ | -2 | -1 | 0 | 1
---|---|---|---|---
up-type | $T^{1}$ | $T^{2}$ | $T^{3}$ | $T^{4}$
|  | $T^{\prime 1}$ | $T^{\prime 2}$ |
down-type | $B^{1}$ | $B^{2}$ | $B^{3}$ |
|  | $B^{\prime 1}$ | $B^{\prime 2}$ |
|  |  | $B^{\prime\prime 1}$ |
Table 2: $U(1)_{\rm{FN}}$ charges for VLQ’s in model 2. Total $6+6=27$ heavy
quarks are needed in this model.

The $U(1)_{\rm{FN}}$ charge assignments are listed in table 2. Clearly, model
2 has a simpler particle content than model 1, although this advantage comes
at the expense of a larger scale separation between $v_{s}$ and $M$. The mass
matrices are,

$\mathcal{M}^{u}=\begin{pmatrix}0&0&0&\langle\widetilde{H}\rangle&0&0&0&0&0\\\
0&0&0&0&\langle\widetilde{H}\rangle&\langle\widetilde{H}\rangle&0&0&0\\\
0&0&\langle\widetilde{H}\rangle&0&0&0&\langle\widetilde{H}\rangle&\langle\widetilde{H}\rangle&0\\\
0&0&0&M&\langle S\rangle&\langle S\rangle&0&0&0\\\ 0&0&\langle
S\rangle&0&M&0&\langle S\rangle&\langle S\rangle&0\\\ 0&0&\langle
S\rangle&0&0&M&\langle S\rangle&\langle S\rangle&0\\\ 0&\langle
S\rangle&0&0&0&0&M&0&\langle S\rangle\\\ 0&\langle
S\rangle&0&0&0&0&0&M&\langle S\rangle\\\ \langle S\rangle&0&0&0&0&0&0&0&M\\\
\end{pmatrix},$ (14a)

$\mathcal{M}^{d}=\begin{pmatrix}0&0&0&\langle H\rangle&0&0&0&0&0\\\
0&0&0&0&\langle H\rangle&\langle H\rangle&0&0&0\\\ 0&0&0&0&0&0&\langle
H\rangle&\langle H\rangle&\langle H\rangle\\\ 0&0&0&M&\langle S\rangle&\langle
S\rangle&0&0&0\\\ 0&0&0&0&M&0&\langle S\rangle&\langle S\rangle&\langle
S\rangle\\\ 0&0&0&0&0&M&\langle S\rangle&\langle S\rangle&\langle S\rangle\\\
\langle S\rangle&\langle S\rangle&\langle S\rangle&0&0&0&M&0&0\\\ \langle
S\rangle&\langle S\rangle&\langle S\rangle&0&0&0&0&M&0\\\ \langle
S\rangle&\langle S\rangle&\langle S\rangle&0&0&0&0&0&M\\\ \end{pmatrix}$ (14b)

###  3.3 Couplings of the UV theory

Each entry of the matrices given in (13) and (14) is accompanied by an
$\mathcal{O}(1)$ coupling $\eta^{u/d}_{ij}$ of the UV theory, that has been
suppressed. We take the bare mass $M$ to be same for all VLQ’s, so that,

$\eta^{u/d}_{ii}=1,~{}\forall~{}i>3.$ (15)

The rest of the yukawa couplings are related to the $\mathcal{O}(1)$ couplings
$y^{u/d}_{ij}$ described in (4) by relations of the form,

$y^{u/d}_{ij}\sim\prod\eta^{u/d}_{kl},$ (16)

where $\eta^{u/d}_{kl}$ are the relevant couplings appearing in the chain
diagrams corresponding to $y^{u/d}_{ij}$. Although these couplings are
$\mathcal{O}(1)$, their exact value is unimportant for our purpose. The
simplest possible choice is to keep all the yukawas unity. Another way is to
choose them from a uniform distribution between say, $[1/3,3]$. Here we take
all yukawas to be equal and parameterize them as,

$\eta^{u/d}_{ij}=y,$ (17)

where $y$, chosen to be a real number, is a scaling factor, which we treat as
a free parameter characterizing the overall yukawa strength. This simplifying
assumption is reasonable, and the factor $y$ will be useful for book-keeping.
In the limit $y=0$, the VLQ’s completely decouple from the scalar sector,
while the perturbativity of yukawas requires that $|y|\leq\sqrt{4\pi}$. We
make a further simplification by choosing $y\geq 0$.

##  4 Effective potential

To study the flavon PT, we need to know the temperature-dependent effective
potential, which is defined in terms of the generating functional of 1PI
connected diagrams. For pedagogic reviews on the finite tmperature effective
potential, see [62, 63, 64].

At zero temperature, the effective potential for a scalar field $\phi$ can be
written in a loop expansion,

$V(\phi)|_{T=0}=V_{0}(\phi)+V_{1}(\phi)+V_{2}(\phi)+\cdots,$ (18)

where the subscript indicates the loop order. The first term on the R.H.S.,
$V_{0}$ is the ‘classical’ or tree-level contribution, while the first quantum
corrections begin at the one-loop level. The temperature dependence on the
other hand, is induced at one-loop level. Hence the one-loop temperature-
dependent effective potential can be written schematically as:

$V(\phi,T)=V_{0}(\phi)+V_{1}(\phi)+V_{1T}(\phi,T)$ (19)

Below, we compute each contribution to the effective potential separately.

###  4.1 Tree-level potential

The Higgs-flavon scalar potential is essentially an extension of the SM by a
complex scalar (CxSM) [65], with a global $U(1)_{\rm{FN}}$ symmetry,

$\displaystyle V_{0}(S,H)$ $\displaystyle=$
$\displaystyle-\mu_{H}^{2}|H|^{2}+\lambda_{H}|H|^{4}+\lambda_{HS}|S|^{2}|H|^{2}$
(20) $\displaystyle-$ $\displaystyle\mu_{S}^{2}|S|^{2}+\lambda_{S}|S|^{4},$

where all parameters are assumed to be real ($CP$ conserving potential). Tree-
level vacuum stability requires,

$\lambda_{H}\geq 0,~{}\lambda_{S}\geq 0.$ (21)

Further, we take $\lambda_{HS}\geq 0.$ Write the fields as:

$H=\begin{pmatrix}G_{+}\\\ \frac{h+iG_{0}}{\sqrt{2}}\\\
\end{pmatrix},~{}~{}S=\frac{1}{\sqrt{2}}s~{}e^{i\rho},$

where $h,~{}s$ are real scalar fields, $G_{0},~{}G_{\pm}$ are SM goldstone
bosons, and $\rho$ is a pseudocalar.

In recent years, it has become popular to associate the pseudoscalar $\rho$
with an Axion-like Particle (ALP), to additionally address the strong $CP$
problem of SM. The pseudoscalar field in this case has been variously called
as the ‘flavorful axion’ [66], the ‘flaxion’ [67] or the ‘axiflavon’ [68]. The
ALP can acquire a tiny mass via the QCD anomaly, making it a pseudo-Nambu
Goldstone (pNG) boson. Alternately, an explicit $U(1)_{\rm{FN}}$ breaking term
such as $V_{0}(S,H)\supset\lambda_{1}(S^{2}+S^{*2})$ can generate a small mass
term for $\rho$, with the mass scale constrained by low energy phenomenology
[40]. The pseudoscalar plays no role in the dynamics of the flavon phase
transition, which is governed by the real scalar field $s$, so we will not
comment on it further.

When written in terms of field components, the potential also has a remnant
$Z_{2}$ symmetry where $s\rightarrow-s$. Spontaneous breaking of such a
discrete symmetry by the VEV of $s$ can lead to formation of domain walls,
which are disallowed by cosmology as they would eventually dominate the energy
density of the universe, and directly contradict the cosmological principle.
An additional $Z_{2}$-breaking term $V_{0}(S,H)\supset\lambda_{2}(S+S^{*})$,
ensures that domain walls disintegrate fast enough to respect cosmological
constraints. The contribution of these explicit symmetry breaking terms is
however small and insignificant for analysing the nature of the flavon PT.

In terms of components, the potential is,

$\displaystyle V_{0}(h,s,G_{0},G_{\pm})$ $\displaystyle=$
$\displaystyle-\frac{\mu^{2}_{H}}{2}\big{(}h^{2}+G_{0}^{2}+2G_{+}G_{-}\big{)}+\frac{\lambda_{H}}{4}\big{(}h^{2}+G_{0}^{2}+2G_{+}G_{-}\big{)}^{2}$
(22) $\displaystyle+$
$\displaystyle\frac{\lambda_{HS}}{4}s^{2}\big{(}h^{2}+G_{0}^{2}+2G_{+}G_{-}\big{)}-\frac{\mu^{2}_{S}}{2}s^{2}+\frac{\lambda_{S}}{4}s^{4}.$

Setting the goldstones to zero,

$V_{0}(h,s)=-\frac{\mu^{2}_{H}}{2}h^{2}+\frac{\lambda_{H}}{4}h^{4}+\frac{\lambda_{HS}}{4}h^{2}s^{2}-\frac{\mu^{2}_{S}}{2}s^{2}+\frac{\lambda_{S}}{4}s^{4}.$
(23)

Suppose at zero temperature the fields get a VEV:

$\langle h\rangle=v_{h},~{}~{}\langle s\rangle=v_{s}.$

We take $\mu^{2}_{H}$, $\lambda_{HS}$, $\lambda_{S}$ and $v_{s}$ as free
parameters. Minimization conditions are,

$\left.\frac{\partial V_{0}}{\partial h}\right|_{(v_{h},v_{s})}=0,$ (24a)
$\left.\frac{\partial V_{0}}{\partial s}\right|_{(v_{h},v_{s})}=0.$ (24b)

These imply,

$\big{(}-\mu^{2}_{H}+\lambda_{H}v_{h}^{2}+\frac{\lambda_{HS}}{2}v_{s}^{2}\big{)}v_{h}=0,$
(25a)
$\big{(}-\mu^{2}_{S}+\lambda_{S}v_{s}^{2}+\frac{\lambda_{HS}}{2}v_{h}^{2}\big{)}v_{s}=0.$
(25b)

Since $v_{s}\neq 0$, the second equation gives,

$\displaystyle\mu_{S}^{2}=\lambda_{S}v_{s}^{2}+\frac{\lambda_{HS}}{2}v_{h}^{2},$
(26)

while in the first equation, we set $v_{h}=0$, so that $\mu_{H}^{2}$ is an
undetermined parameter. The scalar mass matrix is computed from the Hessian of
the potential. Write $(h,s)=(\phi_{1},\phi_{2})$, then,

$m^{2}_{ij}=\left.\frac{\partial^{2}V_{0}}{\partial\phi_{i}\partial\phi_{j}}\right|_{(v_{h},v_{s})}$

Setting $v_{h}=0$,

$\displaystyle\mathcal{M}^{2}_{\rm{scalar}}=\begin{pmatrix}-\mu^{2}_{H}+\frac{\lambda_{HS}}{2}v_{s}^{2}&0\\\
0&2\lambda_{S}v_{s}^{2}\\\ \end{pmatrix}.$ (27)

As physical masses must be real, which gives us an additional constraint,

$\mu^{2}_{H}\leq\frac{\lambda_{HS}}{2}v_{s}^{2}.$ (28)

We make the choice,

$\mu^{2}_{H}=0.01~{}\lambda_{HS}v_{s}^{2},$ (29)

which numerically ensures that if a PT is possible, it occurs along the $s$
direction at one-loop level. This is however not the only possible choice, and
it was verified that changing the multipliying coefficient upto an order of
magnitude doesn’t affect our results significantly. With the above choice, we
can set $h=0$ in the potential in the rest of the analysis.

$V_{0}(s)=-\frac{\mu_{S}^{2}}{2}~{}s^{2}+\frac{\lambda_{S}}{4}~{}s^{4}.$ (30)

The minimization condition (26) now becomes,

$\mu^{2}_{S}=\lambda_{S}v_{s}^{2}.$ (31)

###  4.2 One-loop correction

One-loop correction to the effective potential is given by the Coleman-
Weinberg (CW) formula [69]. In the $\overline{\rm{MS}}$ scheme, in Landau
gauge,

$V_{\rm{CW}}(h,s)=\frac{1}{64\pi^{2}}\sum_{i}(-1)^{f_{i}}n_{i}m^{4}_{i}(h,s)\bigg{[}\log\bigg{(}\frac{m^{2}_{i}(h,s)}{\mu^{2}}\bigg{)}-c_{i}\bigg{]},$
(32)

where the sum runs over all particles that couple to the flavon333We retain
the field $h$ for the sake of completeness, although finally $h$ is to be set
to zero.. These include the SM Higgs ($h$), the SM goldstones
($G_{0},G_{\pm}$), the flavon ($s$), and the SM quarks and VLQ’s ($F$). As the
SM gauge bosons do not directly couple to $s$, they do not contribute to the
one-loop effective potential. In the above equation, $\mu$ is the
renormalization scale, which we set to $\mu^{2}=v_{s}^{2}$. The factor $f_{i}$
is $+1$ $(0)$ for fermions (bosons). $n_{i}$ is the number of degrees of
freedom: (1,1,3,12) for ($h,~{}s,~{}G_{0,{\pm}},~{}F$), and $c_{i}=5/6$ for
gauge bosons and $c_{i}=3/2$ for others.

The field-dependent masses, $m_{i}^{2}(h,s)$, are obtained by expanding the
tree-level Lagrangian around a constant field value $(h,s)$, and reading off
the respective mass terms for the respective species. For the scalars, we get,

$\mathcal{M}^{2}(h,s)=\begin{pmatrix}-\mu^{2}_{H}+3\lambda_{H}h^{2}+\frac{\lambda_{HS}}{2}s^{2}&\lambda_{HS}hs\\\
\lambda_{HS}hs&\lambda_{S}(3s^{2}-v_{s}^{2})+\frac{\lambda_{HS}}{2}h^{2}\\\
\end{pmatrix},$ (33a)
$m^{2}_{G_{0,\pm}}(s)=-\mu^{2}_{H}+\lambda_{H}h^{2}+\frac{\lambda_{HS}}{2}s^{2}.$
(33b)

Note that $m^{2}_{i}(s)<0$ for above scalars near the origin, implying an
imaginary field-dependent mass. This well known problem occurs for scalars
associated with a spontaneously broken symmetry, leading to a complex
potential. The imaginary part of the potential is related to particle decay
[70]. However, the dynamics of PT is governed only by the real part of the
potential and thus we always take the real part of the potential in our
analysis.

For fermions, the field-dependent mass matrices are obtained from the physical
mass matrices $\mathcal{M}^{u/d}$ given in (13) and (14), by the substitution:

$\langle\widetilde{H}\rangle\rightarrow\frac{1}{\sqrt{2}}\begin{pmatrix}h\\\
0\\\ \end{pmatrix},~{}\langle
H\rangle\rightarrow\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\ h\\\
\end{pmatrix},~{}\langle S\rangle\rightarrow\frac{s}{\sqrt{2}}.$

$m_{i}^{2}(h,s)$ are then given by singular values of the field-dependent
$\mathcal{M}^{u/d}$. For matrices of dimension larger than 5, the eigenvalues
need to be evaluated numerically. The lightest three eigenvalues correspond to
the SM quarks, while the rest are the VLQ’s. Contribution of the light quarks
to the potential is negligible and is therefore ignored. In figure 2, we plot
the squared field-dependent masses for the VLQ’s of model 1. In the region
$0\leq s\lesssim M$, with $h=0$, the VLQ masses are well fitted by quadratic
polynomials,

$m^{2}_{i}(s)=M^{2}+a_{i}~{}yMs+b_{i}~{}y^{2}s^{2},$ (34)

Figure 2:  Field-dependent mass squared for the VLQ’s of model 1. Here, we
have taken $y=0.5$ and $v_{s}=10^{5}$ GeV. The cyan coloured vertical line
indicates $s=v_{s}$, where the physical VLQ masses are defined.

where $a_{i}$ and $b_{i}$ are dimensionless fitting parameters, and $y$ is the
scaling factor described in (17). For the largest and smallest VLQ masses, the
parameters were found to be $|a_{i}|,~{}|b_{i}|\sim 1$, while they are smaller
for others.

The $\overline{\rm{MS}}$ renormalization takes care of the UV divergent terms
arising from integrals over momenta running in the loop. However, additional
finite counter-terms can be added to $V_{\rm{CW}}$ so that the VEV and flavon
mass are unchanged from their tree-level values. Write,

$V_{1}(s)=V_{\rm{CW}}(s)+V_{\rm{c.t.}}(s),$ (35a)
$V_{\rm{c.t.}}(s)=-\frac{\delta\mu^{2}_{s}}{2}s^{2}+\frac{\delta\lambda_{s}}{4}s^{4}.$
(35b)

This is the so-called ‘On-shell renormalization’ scheme. From (35a) and (35b),
the finite renormalization conditions are:

$\left.\frac{\partial V_{1}}{\partial s}\right|_{v_{s}}=0,$ (36a)
$\left.\frac{\partial^{2}V_{1}}{\partial s^{2}}\right|_{v_{s}}=0.$ (36b)

These determine the parameters $\delta\mu^{2}_{s}$ and $\delta\lambda_{s}$,

$\delta\mu^{2}_{s}=\frac{1}{2v_{s}}\left(3\left.\frac{\partial
V_{\rm{CW}}}{\partial
s}\right|_{v_{s}}-v_{s}\left.\frac{\partial^{2}V_{\rm{CW}}}{\partial
s^{2}}\right|_{v_{s}}\right),$ (37a)
$\delta\lambda_{s}=\frac{1}{2v_{s}^{3}}\left(\left.\frac{\partial
V_{\rm{CW}}}{\partial
s}\right|_{v_{s}}-v_{s}\left.\frac{\partial^{2}V_{\rm{CW}}}{\partial
s^{2}}\right|_{v_{s}}\right).$ (37b)

The scale-independence of the full effective potential [69] is spoiled by the
appearance of $\mu$ in the CW formula. This is a result of truncating the
potential at one-loop. Approximate scale-independence at one-loop is recovered
when the Logs are $\mathcal{O}(1)$. Hence the renormalization scale needs to
be chosen wisely, depending upon the regime of interest.

In the scenario considered here, there are two scales: $v_{s}$ and the VLQ
mass scale $M$, with the scale separation quantified by the ratio,

$\epsilon\sim\frac{1}{\sqrt{2}}\frac{v_{s}}{M}.$ (38)

For model 1, $\epsilon=0.23$, while for model 2, $\epsilon=0.05$, implying a
larger scale separation for model 2. To study the PT, we are interested in the
behaviour of the potential around $s\sim v_{s}$. For model 1, the largest Log
is $\sim 2.5$, while for model 2, the largest Log is $\sim 5.5$ at
$\mu=v_{s}$. Hence the Logs are $\mathcal{O}(1)$. We treat these Logs as
small, and analyse the one-loop potential at the fixed scale $\mu=v_{s}$, as
the detailed RGE analysis is beyond the scope of this work.

###  4.3 Finite temperature contribution

The one-loop finite temperature correction is essentially the free energy
associated with the scalar field $\phi$ [62],

$V_{1T}(\phi,T)=\sum_{i}n_{i}\frac{T^{4}}{2\pi^{2}}J_{b/f}\bigg{(}\frac{m^{2}_{i}(\phi)}{T^{2}}\bigg{)},$
(39)

where the function $J_{b}$ ($J_{f}$) is defined for bosons (fermions). Write
$x^{2}\equiv\frac{m^{2}_{i}}{T^{2}}$, then $J_{b/f}$ is given by,

$J_{b/f}(x^{2})=\pm\int_{0}^{\infty}dy~{}y^{2}\log\big{(}1\mp
e^{-\sqrt{x^{2}+y^{2}}}\big{)}.$ (40)

These functions have well known expansions for $x^{2}<<1$, and $x^{2}>>1$ .
For small $x^{2}$ [71],

$\displaystyle J_{f}(x^{2},n)=$ $\displaystyle-$
$\displaystyle\frac{7\pi^{4}}{360}+\frac{\pi^{2}}{24}x^{2}+\frac{1}{32}x^{4}\big{(}\log
x^{2}-c_{f}\big{)}$ (41) $\displaystyle-$
$\displaystyle\pi^{2}x^{2}\sum_{l=2}^{n}\left(-\frac{1}{4\pi^{2}}x^{2}\right)^{l}\frac{(2l-3)!!\zeta(2l-1)}{(2l)!!(l+1)}\left(2^{2l-1}-1\right),$
$\displaystyle J_{b}(x^{2},n)=$ $\displaystyle-$
$\displaystyle\frac{\pi^{4}}{45}+\frac{\pi^{2}}{12}x^{2}-\frac{\pi}{6}\big{(}x^{2}\big{)}^{3/2}-\frac{1}{32}x^{4}\big{(}\log
x^{2}-c_{b}\big{)}$ (42) $\displaystyle+$
$\displaystyle\pi^{2}x^{2}\sum_{l=2}^{n}\left(-\frac{1}{4\pi^{2}}x^{2}\right)^{l}\frac{(2l-3)!!\zeta(2l-1)}{(2l)!!(l+1)},$

where, $c_{+}=3/2+2\log\pi-2\gamma_{E}$, and $c_{-}=c_{+}+2\log 4$. Here
$\gamma_{E}$ is the Euler-Mascheroni constant, and $\zeta(x)$ is the
Riemann-$\zeta$ function. $n$ is the order till we want to expand the sum.
Usually taking $n=3,4$ gives good results.

For large $x^{2}$, both fermions and bosons have the same expansion at order
$n$ [71],

$\displaystyle J_{b/f}(x^{2},n)=$
$\displaystyle-\exp\bigg{(}-(x^{2})^{1/2}\bigg{)}\bigg{(}\frac{\pi}{2}(x^{2})^{3/2}\bigg{)}^{1/2}\sum_{l=0}^{n}\frac{1}{2^{l}l!}\frac{\Gamma(5/2+l)}{\Gamma(5/2-l)}(x^{2})^{-l/2},$
(43)

where $\Gamma(x)$ is the Euler Gamma function. For intermediate values of
$x^{2}$, the functions $J_{b/f}(x^{2})$ can be obtained by interpolating
between the large $x^{2}$ and small $x^{2}$ approximation. The finite
temperature effects have been computed in this work using the CosmoTransitions
package [72].

###  4.4 Daisy resummation

The perturbative loop expansion of the effective potential breaks down near
the critical temperature $T_{c}$ due to infrared divergences from the bosonic
Matsubara zero-modes [62]. Leading contribution to the divergences comes from
multiloop ‘daisy diagrams’, which needs to be resummed. We achieve this by
making the substitution: $m^{2}_{i}(s)\rightarrow m^{2}_{i}(s)+\Pi_{i}(T)$
[73], for all the bosons, where $\Pi_{i}(T)$ is the squared thermal mass for
species $i$. Thermal masses are obtained by substituting (33a), (33b) in the
high-$T$ expansion (42) to leading order in $T$. These are:

$\displaystyle\Pi_{h,G_{0,\pm}}(T)$ $\displaystyle=$ $\displaystyle
T^{2}\bigg{(}\frac{g_{1}^{2}}{16}+\frac{3g_{2}^{2}}{16}+\frac{\lambda_{HS}}{12}+\frac{\lambda_{H}}{2}+\frac{y_{t}^{2}}{4}\bigg{)},$
(44) $\displaystyle\Pi_{s}(T)$ $\displaystyle=$ $\displaystyle
T^{2}\bigg{(}\frac{\lambda_{HS}}{6}+\frac{\lambda_{S}}{3}\bigg{)},$ (45)

where $g_{1},~{}g_{2}$ are the usual SM gauge couplings, and $y_{t}$ is the
top yukawa coupling. All SM couplings are taken at the scale $\mu=v_{s}$, and
are found by evolving the SM RGE’s [74].

Now that all the contributions to the effective potential are defined, we
finally have:

$V_{\rm{eff}}(s,T)=V_{0}(s)+V_{\rm{CW}}(s,T)+V_{{\rm{c.t.}}}(s)+V_{1T}(s,T),$
(46)

with the Daisy-corrected field-dependent masses.

##  5 Vacuum stability

Fermions contribute with a negative sign to the CW formula (32), and tend to
destabilize the potential at large field values. In the context of FN, vacuum
stability has been addressed earlier [75]. Similarly, Higgs vacuum stability
of the SM in presence of vector-like fermions has also been discussed [43, 44,
45, 76]. Typically, UV completions of FN contain several ($\mathcal{O}(10)$)
heavy fermions with $\mathcal{O}(1)$ yukawa couplings to the flavon,
destabilizing the effective potential. Here we discuss the vacuum stability
along the flavon direction.

The negative contribution of VLQ’s can be seen from the effective quartic
coupling at one-loop level. Setting the renormalization scale
$\mu^{2}=v_{s}^{2}$ and using (35a),

$\displaystyle\frac{\lambda_{s,\rm{eff}}}{4}$ $\displaystyle=$
$\displaystyle\frac{1}{4!}\left.\frac{\partial^{4}}{\partial
s^{4}}\big{(}V_{0}+V_{1}\big{)}\right|_{s=0}$ (47) $\displaystyle\approx$
$\displaystyle\frac{\lambda_{s}+\delta\lambda_{s}}{4}+\frac{9}{64\pi^{2}}\log\big{(}\lambda_{s}\big{)}+\frac{\lambda^{2}_{HS}}{64\pi^{2}}\log\left(\frac{\lambda_{HS}}{100}\right)$
$\displaystyle+$
$\displaystyle\sum_{i\in{\rm{VLQ}}}\frac{y^{4}}{64\pi^{2}}\bigg{\\{}a_{i}^{4}-12a^{2}_{i}b_{i}-12b_{i}^{2}\log\bigg{(}\frac{M^{2}}{v_{s}^{2}}\bigg{)}\bigg{\\}},$

where the second term on the RHS is the flavon contribution, the third term is
contributed by $h,~{}G_{0,\pm}$, and we have used the fitting polynomials (34)
in the last term for the VLQ contribution. Since the coefficients $a_{i}$ and
$b_{i}$ are all $\mathcal{O}(1)$, the VLQ term is dominated by the factors of
12, and hence the fermionic contribution to the quartic coupling is negative.
Assuming for simplicity that each term in the VLQ sum contributes roughly
equally, the total contribution also increases with the number of VLQ’s.
Lastly, the VLQ contribution scales as $y^{4}$, therefore, increasing the
yukawa factor by a small amount dramatically increases the negative
contribution.

The instability scale, $\Lambda_{\rm{inst}}$ is estimated from the field value
at which the potential falls below the height of the local minimum at
$s=v_{s}$. For $\lambda_{S}\sim\mathcal{O}(0.1)$,
$\lambda_{HS}\sim\mathcal{O}(1)$ and $y=1$, $\Lambda_{\rm{inst}}$ was found
around the VLQ mass scale $M$. Increasing the yukawa scaling factor further
pushes the instability scale lower. Naively, the potential can be stabilized
by increasing the couplings ($\lambda_{S}$, $\lambda_{HS}$), to counteract the
fermionic contribution. However, this procedure is limited by the
perturbativity bounds on the quartic couplings: $\lambda_{S},\lambda_{HS}\leq
4\pi$. For large couplings, the one-loop results also become unreliable, and
any FOPT observed from the one-loop analysis can’t be trusted. As is evident
from (47), the fermion contribution to the quartic coupling is highly
sensitive to the scaling factor $y$. Hence we can push $\Lambda_{\rm{inst}}$
to at least a few times of $M$ by decreasing $y$. We found that by keeping
$y=0.5$, $\Lambda_{\rm{inst}}$ is atleast $~{}3M$ for the smallest value of
the couplings $(\lambda_{S},\lambda_{HS})$ in the parameter region of interest
444Keeping $y=0.5$ for all couplings uniformly can contribute slightly to the
hierarchy of fermion masses as some combinations $y_{ij}$ may be smaller than
$\mathcal{O}(1)$ according to (16), but this effect is tolerable.. Beyond the
instability scale, the potential can be stabilized by adding new heavy bosonic
degrees of freedom which couple to the flavon as illustrated in Appendix B. If
the mass scale of the new bosons is high compared to $M$, we can write down an
Effective Field Theory (EFT), by introducing higher dimensional operators in
the flavon potential in a model-independent way (see for example [43]),

$V_{1}\rightarrow V_{1}+\frac{1}{\Lambda_{b}^{2}}|S|^{6},$ (48)

where $\Lambda_{b}\sim\Lambda_{\rm{inst}}$ is the cut-off scale associated
with the new bosons, and the Wilson coefficient is taken to be unity. For
sufficiently heavy bosons, the thermal contribution is suppressed by the
Boltzmann factor, as can be seen from the small $T$ expansion (43). Hence the
heavy bosonic fields play an insignificant role in determining the nature and
strength of the PT. We assume that vacuum stability at field values beyond
$\Lambda_{\rm{inst}}$ is taken care of by such bosons, without bothering about
any particular construction.

##  6 First order flavon phase transition

###  6.1 Characterizing FOPT

A FOPT proceeds via bubble nucleation, and is characterized by three
parameters related to the latent heat, rate of PT, and the nucleation
temperature, which we review below. As the universe cools and the unbroken
vacuum becomes metastable, tiny bubbles of the stable phase nucleate. In what
follows, we take $T$ as the temperature of the universe at time $t$ after the
Big Bang. The temperature at which the two minima become degenerate is called
the critical temperature $T_{c}$, below which the nucleation rate per unit
volume, $\Gamma(T)$, is given by [77, 78],

$\Gamma(T)=\left(\frac{S_{3}(T)}{2\pi T}\right)^{3/2}T^{4}e^{-S_{3}(T)/T},$
(49)

where $S_{3}(T)$ is the $O(3)$ symmetric Euclidean bounce action, for
transitions at finite temperature. We calculate $S_{3}$ using the Python based
CosmoTransitions [72] package. The tunnelling action $S_{3}$ decreases with
time as the universe cools below $T=T_{c}$, thereby increasing the tunnelling
probability. The nucleation temperature is defined as the temperature $T_{n}$,
at which the tunneling probability per unit volume becomes $\mathcal{O}(1)$,
i.e., at $T=T_{n}$, at least one bubble gets nucleated per Hubble volume. This
happens when,

$\Gamma(T_{n})\approx\big{(}H(T_{n})\big{)}^{4}.$ (50)

In the radiation dominated era, using (49), this implies the condition [58],

$\frac{S_{3}(T)}{T}\simeq-4\ln\left(\frac{T}{m_{\rm{Pl}}}\right),$ (51)

where $m_{\rm{Pl}}$ is the Planck mass.

The strength of FOPT is characterized by the parameter $\alpha$, given by the
ratio of the vacuum energy density and the radiation density of the
surrounding plasma. Let $T_{*}$ be the temperature around which GW production
is most significant. Then,

$\alpha=\left.\frac{\rho_{\rm{vac}}}{\rho_{\rm{rad}}}\right|_{T_{*}},$ (52)

where $\rho_{\rm{vac}}$ is the energy density difference between the
metastable and the stable phases at $T=T_{*}$:
$\rho_{\rm{vac}}=\Delta\bigg{(}V-\frac{\partial V}{\partial T}\bigg{)}$.
$\rho_{\rm{rad}}$ is the radiation density of the thermal bath given by,

$\rho_{\rm{rad}}(T)=\frac{\pi^{2}}{30}g_{*}T^{4}.$ (53)

We choose $T_{*}\approx T_{n}$ in our calculations, since this is where
significant GW production is expected. Around $T_{*}$, $g_{*}=107.75$, with
$106.75$ contributed by the SM fields, and 1 contributed by the pseudoscalar
$\rho$, while the rest of the fields remain non-relativistic.

Finally we need the parameter $\beta$, which characterizes the rate at which
the PT proceeds.

$\beta\equiv-\left.\frac{dS}{dt}\right|_{t=t_{*}}=TH_{*}\left.\frac{dS}{dT}\right|_{T=T_{*}},$
(54)

where, $S=S_{3}/T$, and $H_{*}$ is the Hubble’s constant at $T=T*$. The
negative sign in the first equality is due to the fact that the action
decreases with time.

###  6.2 A case for flavon FOPT

To generate an FOPT, a barrier must be induced between the two degenerate
minima at $T_{c}$. Although the conventional wisdom is that only bosons
contribute to the barrier to generate a SFOPT, it has been shown in the
context of electroweak baryogenesis, that a SFOPT is also possible in purely
fermionic extentions of the SM [79, 80, 81, 82, 83, 84]. It is therefore
interesting to see the effect of the heavy VLQ’s on the strength of FOPT in a
flavon PT. Let us parametrize the finite temperature effective potential,
$V_{\rm{eff}}$ as,

$V_{\rm{eff}}(s,T)=\frac{\mu^{2}(T)}{2}s^{2}+\frac{\lambda_{3}(T)}{3}s^{3}+\frac{\lambda_{4}(T)}{4}s^{4}+\frac{\lambda_{6}(T)}{6}s^{6}.$
(55)

The following two cases may arise:

  * •

Case 1: When the dimension 6 term is negligible, and

$\mu^{2}(T_{c})>0,~{}\lambda_{3}(T_{c})<0,~{}\lambda_{4}(T_{c})>0.$

The cubic term is primarily responsible for inducing the barrier. Such a cubic
term is generated for example by bosons in the high $T$ expansion as in (42),
while fermions do not induce any such contribution. However, both fermions and
bosons have the same expansion in the low $T$ limit (43), and hence can
contribute equally to the formation of a barrier [82].

  * •

Case 2: If the cubic term is small compared to the others, then a barrier can
also be formed when [85],

$\mu^{2}(T_{c})>0,~{}\lambda_{4}(T_{c})<0,~{}\lambda_{6}(T_{c})>0.$

Although the $|S|^{6}$ operator may be generated by integrating out the
stabilizing heavy bosonic fields, we checked that it plays a negligible role
in forming a barrier in our parameter space of interest.

Figure 3:  Shape of the potential of model 1 near the FOPT, for $v_{s}=10^{5}$
GeV, $y=0.5$, $\lambda_{S}=0.105$, $\lambda_{HS}=4.253$.

Clearly the first case is relevant to our scenario, where the effect of a
possible $|S|^{6}$ operator is taken to be small, as discussed in the previous
section. In figure 3, we depict the barrier formed at $T=T_{c}$ (blue), and at
$T=T_{n}$ (orange), for a benchmark point in model 1. The ratio,
$\frac{v_{s}(T_{n})}{T_{n}}\gtrsim 2.7$, indicates a SFOPT. At $T=T_{n}$, the
bosons $\\{h,s,G_{0,\pm}\\}$ have masses $\lesssim T_{n}$, while the VLQ
masses are larger than $T_{n}$. As a result, the bosonic contribution to the
effective potential is governed by the high-$T$ expansion (42), which
incidently justifies the need for taking Daisy resummation into account. On
the other hand, the fermionic contribution is governed by the low-$T$
expansion (43). In fact, this behaviour is observed in the entire region of
parameter space considered in the next section. Therefore, barrier formation
happens primarily by bosons via the cubic term in the high-$T$ expansion, and
the fermionic contribution is comparatively small due to the Boltzmann
suppression factor in (43). In the analysis below, the VLQ’s are observed
reduce the strength of FOPT slightly when compared to the no-VLQ case.

###  6.3 Parameter space

Now we discuss the prospects of achieving a SFOPT in the two FN models (6) and
(9). The parameters considered are, $\\{y,v_{s},\lambda_{S},\lambda_{HS}\\}$.
The VLQ mass scale $M$ is directly related to $v_{s}$ via (38), and is
therefore not an independent parameter.

In what follows, we keep $v_{s}=10^{5}$ GeV, $y\in\\{0,0.5\\}$, and perform a
parameter scan in the $\lambda_{S}-\lambda_{HS}$ plane, with
$\lambda_{S}\in[0,0.3]$, and $\lambda_{HS}\in[2.5,6]$. By restricting
ourselves to this particular region of the plane, we ensure that
perturbativity is obeyed, and that the one-loop calculations are reliable. For
the two models, we choose $y=0.5$, along with the described lower bounds on
$\lambda_{S}$ and $\lambda_{HS}$ so that the instability scale
$\lambda_{\rm{inst}}$ is at least $\gtrsim 3M$. For each point on the plane,
we explore the phase structure of the finite temperature potential using the
‘findallTransitions’ and ‘calcTctrans’ methods of the Python based package
CosmoTransitions [72]. In case of an FOPT, the tunnelling action $S_{3}$ is
computed using the class ‘pathDeformation’, and $T_{n}$ is numerically
obtained using the condition (51). Next, the other two PT parameters, namely
$\alpha$ and $\beta/H_{*}$ are computed. To compute $\beta$, we plot the
action around $T=T_{n}$, and calculate the first derivative using the fourth
order finite difference formula.

Figure 4:  FOPT parameter space when the influence of VLQ’s is switched off.
Colors show PT parameters: $\alpha$ (left), $\beta/H_{*}$ (middle), $T_{n}$
(right). Here, $y=0$, and $v_{s}=10^{5}$ GeV.

First, let us examine the region of FOPT without taking the VLQ’s into
consideration, which corresponds to setting $y=0$ in either model 1 or model
2. Let us call it the ‘reference model’. The region of FOPT is shown in figure
4, where $v_{s}$ is set to $10^{5}$ GeV, as a function of $\alpha$ (left
panel), $\beta/H_{*}$ (middle panel), and $T_{n}$ (right panel). The observed
range of the PT parameters is: $\alpha\in[10^{-3},10^{-1}]$,
$\beta/H_{*}\in[10^{1},10^{6}]$, and $T_{n}\in[10,100]$ TeV. A SFOPT occurs
when $\alpha$ is large, and $\beta/H_{*}$ is small. The region of FOPT is
defined by a sharp boundary on the left hand side of the
$\lambda_{S}-\lambda_{HS}$ plane. In the white region beyond the left
boundary, the metastable minimum at $s=0$, is deep enough, so that the
tunnelling to the $s=v_{s}$ minimum never achieves completion, as the
condition (51) is never met. For a given $\lambda_{HS}$, the boundary occurs
at the lowest value of $\lambda_{S}$ where the FOPT can reach completion, and
this is where it is the strongest. As we move to the right along the
$\lambda_{S}$-axis, the FOPT becomes progressively weaker, as shown by the
plots in figure 4. Beyond the boundary on the right hand side, the FOPT is
either extremely weak, or the PT is second order. Notice that the nucleation
temperature $T_{n}$ increases with increasing couplings $\lambda_{S}$ and
$\lambda_{HS}$.

Figure 5:  PT parameters for model 1 with flavon VEV, $v_{s}=10^{5}$ GeV, and
$y=0.5$. Left panel shows the variation of $\alpha$, with contours
$\alpha=0.012,~{}0.02,~{}0.03,~{}0.04,~{}0.05$. In the middle panel,
$\beta/H_{*}$ is depicted with contours $\beta/H_{*}=3\times
10^{2},~{}10^{3},10^{4},~{}10^{5},~{}10^{6}$. In the right panel, $T_{n}$ is
depicted with contours $T_{n}=4,~{}6,~{}8,~{}9$ $\times 10^{4}$ GeV.

Figure 6:  PT parameters for model 2 with flavon VEV, $v_{s}=10^{5}$ GeV, and
$y=0.5$. $\alpha$ (left), $\beta/H_{*}$ (middle), $T_{n}$ (right). Left panel:
variation of $\alpha$, with contours
$\alpha=0.012,~{}0.02,~{}0.03,~{}0.04,~{}0.05$. Middle panel: $\beta/H_{*}$ is
depicted with contours $\beta/H_{*}=3\times
10^{2},~{}10^{3},10^{4},~{}10^{5},~{}10^{6}$. Right panel: $T_{n}$ is depicted
with contours $T_{n}=4,~{}6,~{}8,~{}9$ $\times 10^{4}$ GeV.

Scans showing the allowed region of FOPT for model 1 are given in figure 5. We
see the same qualitative trends as we move from left to right in the plots of
$\alpha$, $\beta/H_{*}$ and $T_{n}$ as in case of the reference model. The
biggest difference is seen in the left panel depicting $\alpha$, where the
change near the left boundary is more dramatic compared to the reference
model. The effect of the VLQ’s of model 1 is therefore to reduce the region of
SFOPT slightly. The color gradient in the plots of $\beta/H_{*}$ and $T_{n}$
is similarly, more dramatic as compared to that of the reference model. To
note the differences between model 1 and model 2, we overlay contours for a
few chosen values of $\alpha$, $\beta/H_{*}$ and $T_{n}$.

In figure 6, we show the parameter scans for model 2, with the same choice of
contours as model 1. In the left panel depicting $\alpha$, we see that the
contrours for $\alpha=0.012,~{}0.02$ appear more spread out when compared to
the respective contours for model 1, indicating a slower fall in the $\alpha$
values in model 2 as we increase $\lambda_{S}$. Close to the left boundary
however, the contours of the two models are comparable. We find that the plots
resemble those of the reference plots more closely as compared to model 1.
This is because of the larger scale separation between $v_{s}$ and the VLQ
mass scale $M$ for model 2, compared to model 1, due to which the VLQ
contribution is more suppressed in model 2. As a result, the influence of
VLQ’s is less pronounced in model 2 as compared to model 1. In the next
section, we describe the GW signature arising from a SFOPT.

##  7 Gravitational wave background

The GW spectrum is defined in terms of the GW energy density $\rho_{\rm{GW}}$,
as [28],

$\Omega_{\rm{GW}}(f)\equiv\frac{1}{\rho_{c}}\frac{d\rho_{\rm{GW}}}{d\ln f},$
(56)

where $f$ is the frequency and $\rho_{c}$ is the critical energy density of
the universe, given by,

$\rho_{c}=\frac{3H_{0}^{2}}{8\pi G}.$ (57)

Here, $H_{0}=100h~{}{\rm{km/s/Mpc}}$, is the Hubble constant today, with the
current value of $h=0.675\pm 0.005$, and $G$ is the Newton’s gravitational
constant.

Given a strong FOPT, bubbles of the stable phase larger than a critical size
expand rapidly. If there is sufficient friction in the plasma, the bubble
walls may end up reaching a terminal velocity $v_{w}$. GW’s are produced when
the bubbles collide and coalesce with each other. The latent heat of the
colliding bubbles gets distributed into GW’s via three main processes: bubble
wall collisions ($\Omega_{\rm{col}}$), sound waves produced in the thermal
plasma ($\Omega_{\rm{sw}}$), and the resulting MHD turbulence
($\Omega_{\rm{turb}}$). The contributions can be added in the linear
approximation:

$h^{2}\Omega_{\rm{GW}}\simeq
h^{2}\Omega_{\rm{col}}+h^{2}\Omega_{\rm{sw}}+h^{2}\Omega_{\rm{turb}}.$ (58)

In this paper, we have considered the spontaneous breaking of a global
$U(1)_{\rm{FN}}$ symmetry, where no gauge bosons acquire a mass when they
cross the bubble walls from the symmetric phase to the broken phase. Also,
there are no other light particles in the symmetric phase which become heavy
in the broken phase. As a result, the friction exerted by the plasma on the
bubble walls is expected to be negligible, which means that the walls can
accelerate without any bounds. The runaway scenario [86, 87] is therefore
relevant, and we take $v_{w}=1$. The contribution from bubble walls is
expected to be the most significant while the effect of sound waves and MHD
turbulence is subdominant. For the sake of completeness though, we list all
three contributions to $h^{2}\Omega_{\rm{GW}}$ [28] in the runaway scenario ,

$\displaystyle h^{2}\Omega_{\rm{col}}(f)$ $\displaystyle=$ $\displaystyle
1.67\times
10^{-5}\left(\frac{H_{*}}{\beta}\right)^{2}\left(\frac{\kappa_{\phi}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{1/3}\left(\frac{0.11v_{w}^{3}}{0.42+v_{w}^{2}}\right)S_{\rm{col}}(f),$
(59) $\displaystyle h^{2}\Omega_{\rm{sw}}(f)$ $\displaystyle=$ $\displaystyle
2.65\times
10^{-6}\left(\frac{H_{*}}{\beta}\right)^{2}\left(\frac{\kappa_{v}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{1/3}v_{w}~{}S_{\rm{sw}}(f),$
(60) $\displaystyle h^{2}\Omega_{\rm{turb}}(f)$ $\displaystyle=$
$\displaystyle 3.35\times
10^{-4}\left(\frac{H_{*}}{\beta}\right)^{2}\left(\frac{\kappa_{\rm{turb}}\alpha}{1+\alpha}\right)^{3/2}\left(\frac{100}{g_{*}}\right)^{1/3}v_{w}~{}S_{\rm{turb}}(f).$
(61)

The spectral shape functions, $S_{\rm{col}}$, $S_{\rm{sw}}$ and
$S_{\rm{turb}}$ determine the power law behaviour of each contribution at low
and high frequencies. These are,

$\displaystyle S_{\rm{col}}(f)$ $\displaystyle=$
$\displaystyle\left(\frac{f}{f_{\rm{col}}}\right)^{2.8}\frac{3.8}{1+2.8(f/f_{\rm{col}})^{3.8}},$
(62) $\displaystyle S_{\rm{sw}}(f)$ $\displaystyle=$
$\displaystyle\left(\frac{f}{f_{\rm{sw}}}\right)^{3}\left(\frac{7}{4+3(f/f_{\rm{sw}})^{2}}\right)^{7/2},$
(63) $\displaystyle S_{\rm{turb}}(f)$ $\displaystyle=$
$\displaystyle\left(\frac{f}{f_{\rm{turb}}}\right)^{3}\frac{1}{[1+(f/f_{\rm{turb}})]^{11/3}(1+8\pi
f/h_{*})}.$ (64)

Here, $h_{*}$ is the Hubble rate at $T=T_{*}$,

$h_{*}=1.65\times
10^{-7}~{}{\rm{Hz}}\left(\frac{T_{*}}{100~{}{\rm{GeV}}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}.$
(65)

Taking into account the expansion of the universe from $T=T_{*}$ to present
day, the red-shifted peak frequencies are,

$\displaystyle f_{\rm{col}}$ $\displaystyle=$ $\displaystyle 1.65\times
10^{-5}~{}{\rm{Hz}}\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100~{}{\rm{GeV}}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}\left(\frac{0.62}{1.8-0.1v_{w}+v_{w}^{2}}\right),$
(66) $\displaystyle f_{\rm{sw}}$ $\displaystyle=$ $\displaystyle 1.9\times
10^{-5}~{}{\rm{Hz}}~{}\frac{1}{v_{w}}\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100~{}{\rm{GeV}}}\right)\left(\frac{g_{*}}{100}\right)^{1/6},$
(67) $\displaystyle f_{\rm{turb}}$ $\displaystyle=$ $\displaystyle 2.7\times
10^{-5}~{}{\rm{Hz}}~{}\frac{1}{v_{w}}\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100~{}{\rm{GeV}}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}.$
(68)

We also need the efficiency factors: $\kappa_{\phi}$, $\kappa_{\rm{sw}}$,
$\kappa_{\rm{turb}}$ representing the fraction of the latent heat which is
shared by the three different processes. In the runaway scenario, these are,

$\kappa_{\phi}=\frac{\alpha-\alpha_{\infty}}{\alpha},~{}~{}\kappa_{v}=\frac{\alpha_{\infty}}{\alpha}\kappa_{\infty},~{}~{}\kappa_{\rm{turb}}=\epsilon_{\rm{MHD}}\kappa_{v}.$
(69)

Here, $\epsilon_{\rm{MHD}}$ is the turbulent fraction of bulk motion in
plasma, which is at most $5\%$-$10\%$ [88]. Also, $\kappa_{\infty}$ is given
by,

$\displaystyle\kappa_{\infty}$ $\displaystyle\equiv$
$\displaystyle\frac{\alpha_{\infty}}{0.73+0.083\sqrt{\alpha_{\infty}}+\alpha_{\infty}},$
(70) $\displaystyle\alpha_{\infty}$ $\displaystyle\simeq$
$\displaystyle\frac{5}{4\pi^{2}}\frac{\sum_{i}c_{i}\Delta
m_{i}^{2}}{g_{*}T_{*}^{2}},$ (71)

where the sum over $i$ runs over all particles that are light in the symmetric
phase and acquire mass in the broken phase, with $\Delta m_{i}^{2}$ as the
mass squared difference in the two phases. Also, $c_{i}=n_{i}(n_{i}/2)$ for
bosons (fermions) with $n_{i}$ the number of degrees of freedom of the
particle. Since in our case there are no such light particles which become
heavy as they cross the bubble wall, we expect $\alpha_{\infty}$ to be
negligible. As a result $\alpha>>\alpha_{\infty}$, and only the contribution
from bubble collisions is significant,

$h^{2}\Omega_{\rm{GW}}\simeq h^{2}\Omega_{\rm{col}}.$ (72)

###  7.1 Detection prospects

To get a detectable GW signal resulting from FOPT, we need large $\alpha$, and
small $\beta/H_{*}$, as is evident from the dependence of expressions (59),
(60) and (61) on these quantities. In the parameter scans considered in
section 6, we saw that the strongest FOPT occurs near the left boundary as
shown in figure 4, figure 5, and figure 6.

Figure 7:  GW spectra for model 1 and model 2 for $\lambda_{s}\in[0.1,0.18]$.
Here $\lambda_{HS}=4.253$, $v_{s}=10^{5}$ GeV, and $y=0.5$.

Figure 8:  GW spectra for model 1 and model 2 for
$v_{s}=10^{4,5,6,7}~{}{\rm{GeV}}$. Here $\lambda_{s}=0.105$,
$\lambda_{HS}=4.253$, and $y=0.5$. The peak frequency is seen to shift from
left to right as $v_{s}$ increases.

An example of GW spectra highlighting the dependence on model parameters, is
shown in figure 7. For model 1 and model 2, the dashed curves show the GW
spectrum, when $\lambda_{HS}=4.253$, and $\lambda_{S}$ is varied from
$\lambda_{S}=0.1$ to $\lambda_{S}$ =0.18. The other parameters are kept fixed
at $v_{s}=10^{5}$ GeV, and $y=0.5$. Solid curves are the power law integrated
sensitivity curves (PLISC) [89] for the GW observatories LISA [3], BBO [4],
DECIGO [5], ET [7] and CE [6]. The strength of the GW spectrum decreases as
$\lambda_{S}$ increases. This corresponds to moving horizontally from left to
right, starting close to the boundary in the $\lambda_{S}$-$\lambda_{HS}$
plane in figures 5 and 6. As the GW strength decreases, it eventually goes
below the noise floor of DECIGO and BBO, becoming undetectable.

In figure 8, we show the effect of varying the scale $v_{s}$ on the GW
spectrum for a benchmark value of $(\lambda_{S},\lambda_{HS})$. Taking
$v_{s}=10^{4},~{}10^{5},~{}10^{6},~{}10^{7}~{}{\rm{GeV}}$, we observe that the
peak frequency shifts from left to right as $v_{s}$ increases. This can be
explained from (66),(67) and (68), where the peak frequency is proportional to
$T_{*}\sim v_{s}$. The height of the peak amplitude is seen to increase with
$v_{s}$ due to the temperature dependence of the RHS of (51), which leads to a
larger amount of supercooling at higher $T$. Clearly, different detectors are
suitable for different scales $v_{s}$. For $v_{s}=10^{4,5}$ GeV, the chosen
benchmark point gives a strong GW signal at DECIGO and BBO for both model 1
and model 2. While for $v_{s}=10^{6}$ GeV, the signal would also be detectable
at CE, in addition to DECIGO and BBO. For $v_{s}=10^{7}$ GeV, the signal can
be detected at BBO, DECIGO, CE and ET. For the entire viable region of FOPT,
the GW signal is not observable at LISA, since the peak amplitude is almost
always found to be smaller than $10^{-12}$.

Figure 9:  Regions of GW detection for model 1 and model 2, with
$v_{s}=10^{5}$ GeV, and $y=0.5$. The region coloured in purple overlaps the
yellow region, since BBO is more sensitive compared to DECIGO.

Figure 10:  Regions of GW detection for model 1 and model 2, with
$v_{s}=10^{7}$ GeV, and $y=0.5$. The region coloured in green overlaps the
peach region, since CE is more sensitive compared to ET.

Next, in figure 9, we showcase the detectable regions of the parameter space
for the model 1 (left panel) and model 2 (right panel), when $v_{s}=10^{5}$
GeV, and $y=0.5$. The coloured region corresponds to the region sensitive to
DECIGO (purple), and BBO (yellow). The two regions overlap each other; the
region corresponding to BBO is broader than that of DECIGO, since the former
is more sensitive. In the blue region, the GW signal is too weak to be
detected. In figure 10, we show the detectable parameter space for both models
with $v_{s}=10^{7}$ GeV, and $y=0.5$. Now the parameter space is accessible to
ET and CE in addition to DECIGO and BBO. For convenience we show only regions
for ET (green) and CE (peach), as the regions for BBO and DECIGO are much
smaller. The detectable region is narrower compared to the region at $10^{5}$
GeV. The best detection prospects for $v_{s}=10^{7}$ GeV are at CE.

Previously in figure 5 and figure 6, it was found that the constant $\alpha$
contours for chosen values lie closer to the boundary for model 1 as compared
to model 2. However it is interesting to note that apart from minor
qualitative differences, the GW signal strength is comparable for the two
models in the detectable region, both at $v_{s}=10^{5}$ GeV, and
$v_{s}=10^{7}$ GeV. This is because the regions of detectability require
$\alpha\gtrsim 0.03$, where the contours in the two models are similar. This
implies that a detectable GW background cannot discrimitate between model 1
and model 2. The GW signal is therefore not sensitive to the specific
$U(1)_{\rm{FN}}$ charge assignment, as long as there is only one flavon which
breaks the symmetry.

##  8 Conclusions

In this paper, we have explored the possibility of probing FN models with GW
detectors, if the flavon undergoes a SFOPT. To do so, we have constructed two
minimal, non-supersymmetric UV completions of the FN mechanism involving a
single flavon. To avoid any possible complications arising from gauging it, we
have chosen the family symmetry group as global $U(1)$. However, the two
models considered here are UV complete only in the sense that they generate
the effective FN operators which give rise to the hierarchy of quark yukawas
at low scale. We find that it is essential to add heavy bosons couple to the
flavon, to counter the destabilizing effect of adding several fermions to the
theory. We have chosen a region of parameter space where the effect of the
heavy bosons on the nature and strength of FOPT is insignificant, to make
predictions independent of any specific realization of heavy bosons.

To avoid stringent constraints from collider experiments, we have focussed on
the intermediate energy scale, $10^{4}$-$10^{7}$ GeV for the scale of the FN
mechanism. We find that in the $\lambda_{S}$-$\lambda_{HS}$ plane, a strong
FOPT is observed near the boundary separating the FOPT region from the region
where the PT never completes. The region for observing a GW signal therfore
lies within a thin strip close to this boundary. Detection prospects are the
best when $v_{s}\sim 10^{4}$-$10^{5}$ GeV, in which case the peak of the GW
signal lies in the sensitivity range of BBO and DECIGO. Higher scales of
$v_{s}\sim 10^{6}$-$10^{7}$ GeV can be probed by BBO, DECIGO, as well as CE
and ET.

Detecting a stochastic GW background can become a powerful tool to probe the
flavor symmetry-breaking scale in the coming years, along with complementary
collider signatures. For the two models considered, the GW signatures are of
comparable strength when the same set of parameters is used. Therefore we
conclude that it is not possible to discriminate between them, based solely on
the GW signature. We take the liberty to extrapolate this finding to all
models which involve a single flavon, when the family symmetry group is
$U(1)_{\rm{FN}}$. As such, models involving more than one flavon may be
discriminated by GW’s since they may feature multiple peaks, owing to
different scales of the VEV’s, but this possibility remains to be seen. It
would be interesting to see if the GW signature can throw some light on the FN
symmetry group.

###### Acknowledgements.

D.R. would like to thank Subhendu Rakshit, Siddhartha Karmakar, Dipankar Das,
and Najimuddin Khan for useful discussions. This work is supported by the
Department of Science and Technology, India via SERB grant EMR/2014/001177 and
DST-DAAD grant INT/FRG/DAAD/P-22/2018. D.R. is supported by the Govt. of India
UGC-SRF fellowship.

##  Appendix A A two-flavon model

Consider the horizontal group: $G_{\rm{FN}}=U(1)_{\rm{H1}}\times
U(1)_{\rm{H2}}$. The flavor symmery is then broken by two flavons
$S_{1},~{}S_{2}$, one for each $U(1)$ group, which get a VEV. A possible two-
flavon model is given below [36]:

  * •

Model 3:
$\epsilon_{1}\sim\lambda^{2}=0.04,~{}\epsilon_{2}\sim\lambda^{3}=0.008$,

$Q_{\rm{FN}}(\overline{Q}_{L})=((0,1),(1,0),(0,0)),$
$Q_{\rm{FN}}(u_{R})=((0,1),(-1,1),(0,0)),$
$Q_{\rm{FN}}(d_{R})=((0,1),(1,0),(1,0))$

The yukawa matrices are:

$\displaystyle{\bf{Y^{u}}}\sim\begin{pmatrix}\epsilon_{2}^{2}&\epsilon_{1}^{-1}\epsilon_{2}^{2}&\epsilon_{2}\\\
\epsilon_{1}\epsilon_{2}&\epsilon_{2}&\epsilon_{1}\\\
\epsilon_{2}&\epsilon_{1}^{-1}\epsilon_{2}&1\end{pmatrix},~{}~{}~{}{\bf{Y^{d}}}\sim\begin{pmatrix}\epsilon_{2}^{2}&\epsilon_{1}\epsilon_{2}&\epsilon_{1}\epsilon_{2}\\\
\epsilon_{1}\epsilon_{2}&\epsilon_{1}^{2}&\epsilon_{1}^{2}\\\
\epsilon_{2}&\epsilon_{1}&\epsilon_{1}\end{pmatrix}.$ (73)

The determinants are:

$\det{\bf{Y^{u}}}\sim\epsilon_{2}^{3},~{}~{}~{}\det{\bf{Y^{d}}}\sim\epsilon_{2}^{2}\epsilon_{1}^{3}$
(74)

In supersymmetric models, negative powers of $\epsilon$ are prohibited due to
holomorphy; this gives rise to texture zeros (see [36]). In this work, we have
only considered non-supersymmetric realizations of FN.

##  Appendix B Stabilizing the potential with heavy bosons

Here we illustrate a scenario under which the effective potential can be
stabilized by heavy bosons. Consider $N$ bosons at mass scale $M_{b}$, with an
unbroken $O(N)$ symmetry:

$\displaystyle\mathcal{L}$ $\displaystyle\supset$
$\displaystyle\frac{1}{2}\big{(}\partial^{\mu}\Phi\big{)}^{\dagger}\big{(}\partial_{\mu}\Phi\big{)}-V(\Phi,S)$
$\displaystyle V(\Phi,S)$ $\displaystyle=$
$\displaystyle\frac{1}{2}M_{b}^{2}\big{(}\Phi^{\dagger}\Phi\big{)}+\frac{1}{4}\lambda_{b}\big{(}\Phi^{\dagger}\Phi\big{)}^{2}+\frac{1}{2}\lambda_{bS}\big{(}\Phi^{\dagger}\Phi\big{)}|S|^{2},$
(75)

where $\Phi$ is a column vector of $N$ real scalar fields, and $\lambda_{b}$
and $\lambda_{bS}$ are the new couplings. Including the Coleman-Weinberg
correction from these bosons modifies $V_{\rm{CW}}$ as,

$V_{\rm{CW}}=\frac{1}{64\pi^{2}}\sum_{i}(-1)^{f_{i}}n_{i}m^{4}_{i}\bigg{[}\log\bigg{(}\frac{m^{2}_{i}}{\mu^{2}}\bigg{)}-c_{i}\bigg{]}+\frac{1}{64\pi^{2}}\sum_{b}m^{4}_{i}\bigg{[}\log\bigg{(}\frac{m^{2}_{b}}{\mu^{2}}\bigg{)}-c_{b}\bigg{]},$
(76)

where the sum over $i$ includes SM fields, VLQ’s and the flavon, while sum
over $b$ includes the new heavy bosons whose field dependent mass $m^{2}_{b}$
is given by,

$m^{2}_{b}(s)=M_{b}^{2}+\frac{1}{2}\lambda_{bS}s^{2}.$ (77)

We take $\mu^{2}=v_{s}^{2}$, implying $m_{b}^{2}>>\mu^{2}$, and hence the
heavy bosons can be integrated out, with their effect captured by higher
dimensional operators (see for example [43]). Keeping the leading order EFT
operator, (78) becomes,

$V_{\rm{CW}}=\frac{1}{64\pi^{2}}\sum_{i}(-1)^{f_{i}}n_{i}m^{4}_{i}\bigg{[}\log\bigg{(}\frac{m^{2}_{i}}{\mu^{2}}\bigg{)}-c_{i}\bigg{]}+\frac{c_{6}}{\Lambda_{b}^{2}}s^{6},$
(78)

with $c_{6}$ as the Wilson coefficient.

## References

  * [1] B. P. Abbott et al.  Observation of gravitational waves from a binary black hole merger.  Physical Review Letters, 116(6), feb 2016.
  * [2] Gregory M. Harry.  Advanced LIGO: The next generation of gravitational wave detectors.  Class. Quant. Grav., 27:084006, 2010.
  * [3] Amaro-Seoane et al.  Laser interferometer space antenna, 2017.
  * [4] Vincent Corbin and Neil J Cornish.  Detecting the cosmic gravitational wave background with the big bang observer.  Classical and Quantum Gravity, 23(7):2435–2446, mar 2006.
  * [5] Mitsuru Musha.  Space gravitational wave detector DECIGO/pre-DECIGO.  Proc. SPIE Int. Soc. Opt. Eng., 10562:105623T, 2017.
  * [6] B P Abbott et al.  Exploring the sensitivity of next generation gravitational wave detectors.  Classical and Quantum Gravity, 34(4):044001, jan 2017.
  * [7] M. Punturo et al.  The Einstein Telescope: A third-generation gravitational wave observatory.  Class. Quant. Grav., 27:194002, 2010.
  * [8] Zaven Arzoumanian et al.  The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background.  Astrophys. J. Lett., 905(2):L34, 2020.
  * [9] R. N. Manchester.  The International Pulsar Timing Array.  Class. Quant. Grav., 30:224010, 2013.
  * [10] Sarah Burke-Spolaor et al.  The Astrophysics of Nanohertz Gravitational Waves.  Astron. Astrophys. Rev., 27(1):5, 2019.
  * [11] G. Hobbs et al.  The international pulsar timing array project: using pulsars as a gravitational wave detector.  Class. Quant. Grav., 27:084013, 2010.
  * [12] L. Lentati et al.  European Pulsar Timing Array Limits On An Isotropic Stochastic Gravitational-Wave Background.  Mon. Not. Roy. Astron. Soc., 453(3):2576–2598, 2015.
  * [13] R. N. Manchester et al.  The Parkes Pulsar Timing Array Project.  Publ. Astron. Soc. Austral., 30:17, 2013.
  * [14] F Acernese et al.  Advanced virgo: a second-generation interferometric gravitational wave detector.  Classical and Quantum Gravity, 32(2):024001, dec 2014.
  * [15] B. P. Abbott et al.  GWTC-1: A gravitational-wave transient catalog of compact binary mergers observed by LIGO and virgo during the first and second observing runs.  Physical Review X, 9(3), sep 2019.
  * [16] B. P. Abbott et al.  GW190425: Observation of a compact binary coalescence with total mass $\sim~{}3.4m$.  The Astrophysical Journal Letters, 892(1):L3, mar 2020.
  * [17] Sabino Matarrese, Silvia Mollerach, and Marco Bruni.  Second order perturbations of the Einstein-de Sitter universe.  Phys. Rev. D, 58:043504, 1998.
  * [18] Silvia Mollerach, Diego Harari, and Sabino Matarrese.  CMB polarization from secondary vector and tensor modes.  Phys. Rev. D, 69:063002, 2004.
  * [19] Daniel Baumann, Paul J. Steinhardt, Keitaro Takahashi, and Kiyotomo Ichiki.  Gravitational Wave Spectrum Induced by Primordial Scalar Perturbations.  Phys. Rev. D, 76:084019, 2007.
  * [20] Kazunori Kohri and Takahiro Terada.  Semianalytic calculation of gravitational wave spectrum nonlinearly induced from primordial curvature perturbations.  Phys. Rev. D, 97(12):123532, 2018.
  * [21] Francesco D’Eramo and Kai Schmitz.  Imprint of a scalar era on the primordial spectrum of gravitational waves.  Phys. Rev. Research., 1:013010, 2019.
  * [22] Alexander Vilenkin.  Cosmic Strings and Domain Walls.  Phys. Rept., 121:263–315, 1985.
  * [23] R. R. Caldwell and Bruce Allen.  Cosmological constraints on cosmic string gravitational radiation.  Phys. Rev. D, 45:3447–3468, 1992.
  * [24] M. B. Hindmarsh and T. W. B. Kibble.  Cosmic strings.  Rept. Prog. Phys., 58:477–562, 1995.
  * [25] John Preskill, Sandip P. Trivedi, Frank Wilczek, and Mark B. Wise.  Cosmology and broken discrete symmetry.  Nucl. Phys. B, 363:207–220, 1991.
  * [26] Marcelo Gleiser and Ronald Roberts.  Gravitational waves from collapsing vacuum domains.  Phys. Rev. Lett., 81:5497–5500, 1998.
  * [27] Takashi Hiramatsu, Masahiro Kawasaki, and Ken’ichi Saikawa.  On the estimation of gravitational wave spectrum from cosmic domain walls.  JCAP, 02:031, 2014.
  * [28] Chiara Caprini et al.  Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions.  JCAP, 04:001, 2016.
  * [29] Chiara Caprini and Daniel G. Figueroa.  Cosmological Backgrounds of Gravitational Waves.  Class. Quant. Grav., 35(16):163001, 2018.
  * [30] Georges Aad et al.  Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC.  Phys. Lett. B, 716:1–29, 2012.
  * [31] Serguei Chatrchyan et al.  Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC.  Phys. Lett. B, 716:30–61, 2012.
  * [32] K. S. Babu.  TASI Lectures on Flavor Physics.  In Theoretical Advanced Study Institute in Elementary Particle Physics: The Dawn of the LHC Era, pages 49–123, 2010.
  * [33] Yuval Grossman and Philip Tanedo.  Just a Taste: Lectures on Flavor Physics.  In Theoretical Advanced Study Institute in Elementary Particle Physics: Anticipating the Next Discoveries in Particle Physics, pages 109–295, 2018.
  * [34] Y. Nir.  Flavour physics and CP violation.  CERN Yellow Rep. School Proc., 5:79–128, 2020.
  * [35] C. D. Froggatt and Holger Bech Nielsen.  Hierarchy of Quark Masses, Cabibbo Angles and CP Violation.  Nucl. Phys. B, 147:277–298, 1979.
  * [36] Miriam Leurer, Yosef Nir, and Nathan Seiberg.  Mass matrix models.  Nuclear Physics B, 398(2):319–342, Jun 1993.
  * [37] Miriam Leurer, Yosef Nir, and Nathan Seiberg.  Mass matrix models: the sequel.  Nuclear Physics B, 420(3):468–504, Jun 1994.
  * [38] Lorenzo Calibbi, Zygmunt Lalak, Stefan Pokorski, and Robert Ziegler.  The Messenger Sector of SUSY Flavour Models and Radiative Breaking of Flavour Universality.  JHEP, 06:018, 2012.
  * [39] Leon M. G. de la Vega, Newton Nath, Stefan Nellen, and Eduardo Peinado.  Flavored axion in the UV-complete Froggatt–Nielsen models.  Eur. Phys. J. C, 81(7):608, 2021.
  * [40] Martin Bauer, Torben Schell, and Tilman Plehn.  Hunting the flavon.  Physical Review D, 94(5), Sep 2016.
  * [41] Sebastian A. R. Ellis, Rohini M. Godbole, Shrihari Gopalakrishna, and James D. Wells.  Survey of vector-like fermion extensions of the Standard Model and their phenomenological implications.  JHEP, 09:130, 2014.
  * [42] Mathieu Buchkremer and Alexander Schmidt.  Long-lived heavy quarks : a review.  Adv. High Energy Phys., 2013:690254, 2013.
  * [43] Kfir Blum, Raffaele Tito D’Agnolo, and JiJi Fan.  Vacuum stability bounds on Higgs coupling deviations in the absence of new bosons.  JHEP, 03:166, 2015.
  * [44] Shrihari Gopalakrishna and Arunprasath Velusamy.  Higgs vacuum stability with vectorlike fermions.  Phys. Rev. D, 99(11):115020, 2019.
  * [45] Debasish Borah, Rishav Roshan, and Arunansu Sil.  Sub-TeV singlet scalar dark matter and electroweak vacuum stability with vectorlike fermions.  Phys. Rev. D, 102(7):075034, 2020.
  * [46] The International Linear Collider Technical Design Report - Volume 1: Executive Summary.  6 2013.
  * [47] The International Linear Collider Technical Design Report - Volume 2: Physics.  6 2013.
  * [48] A. Abada et al.  FCC Physics Opportunities: Future Circular Collider Conceptual Design Report Volume 1.  Eur. Phys. J. C, 79(6):474, 2019.
  * [49] A. Abada et al.  FCC-ee: The Lepton Collider: Future Circular Collider Conceptual Design Report Volume 2.  Eur. Phys. J. ST, 228(2):261–623, 2019.
  * [50] A. Abada et al.  FCC-hh: The Hadron Collider: Future Circular Collider Conceptual Design Report Volume 3.  Eur. Phys. J. ST, 228(4):755–1107, 2019.
  * [51] P. S. Bhupal Dev and A. Mazumdar.  Probing the Scale of New Physics by Advanced LIGO/VIRGO.  Phys. Rev. D, 93(10):104001, 2016.
  * [52] Andrea Addazi, Antonino Marcianò, and Roman Pasechnik.  Probing Trans-electroweak First Order Phase Transitions from Gravitational Waves.  MDPI Physics, 1(1):92–102, 2019.
  * [53] Djuna Croon, Tomás E. Gonzalo, and Graham White.  Gravitational Waves from a Pati-Salam Phase Transition.  JHEP, 02:083, 2019.
  * [54] Iason Baldes and Géraldine Servant.  High scale electroweak phase transition: baryogenesis \& symmetry non-restoration.  JHEP, 10:053, 2018.
  * [55] P. S. Bhupal Dev, Francesc Ferrer, Yiyang Zhang, and Yongchao Zhang.  Gravitational Waves from First-Order Phase Transition in a Simple Axion-Like Particle Model.  JCAP, 11:006, 2019.
  * [56] Admir Greljo, Toby Opferkuch, and Ben A. Stefanek.  Gravitational Imprints of Flavor Hierarchies.  Phys. Rev. Lett., 124(17):171802, 2020.
  * [57] Benedict Von Harling, Alex Pomarol, Oriol Pujolàs, and Fabrizio Rompineve.  Peccei-Quinn Phase Transition at LIGO.  JHEP, 04:195, 2020.
  * [58] Wei-Chih Huang, Francesco Sannino, and Zhi-Wei Wang.  Gravitational Waves from Pati-Salam Dynamics.  Phys. Rev. D, 102(9):095025, 2020.
  * [59] Chong-Sun Chu and Asuka Ito.  Gravitational waves in metastable supersymmetry breaking.  Phys. Rev. D, 105(12):123538, 2022.
  * [60] Pierre Binetruy, Stephane Lavignac, and Pierre Ramond.  Yukawa textures with an anomalous horizontal Abelian symmetry.  Nucl. Phys. B, 477:353–377, 1996.
  * [61] Micha Berkooz, Yosef Nir, and Tomer Volansky.  Baryogenesis from the Kobayashi-Maskawa phase.  Phys. Rev. Lett., 93:051301, 2004.
  * [62] M. E. Carrington.  The Effective potential at finite temperature in the Standard Model.  Phys. Rev. D, 45:2933–2944, 1992.
  * [63] Mariano Quiros.  Finite temperature field theory and phase transitions.  In ICTP Summer School in High-Energy Physics and Cosmology, pages 187–259, 1 1999.
  * [64] Mikko Laine and Aleksi Vuorinen.  Basics of Thermal Field Theory, volume 925.  Springer, 2016.
  * [65] Matthew Gonderinger, Hyungjun Lim, and Michael J. Ramsey-Musolf.  Complex Scalar Singlet Dark Matter: Vacuum Stability and Phenomenology.  Phys. Rev. D, 86:043511, 2012.
  * [66] Frank Wilczek.  Axions and Family Symmetry Breaking.  Phys. Rev. Lett., 49:1549–1552, 1982.
  * [67] Yohei Ema, Koichi Hamaguchi, Takeo Moroi, and Kazunori Nakayama.  Flaxion: a minimal extension to solve puzzles in the standard model.  JHEP, 01:096, 2017.
  * [68] Lorenzo Calibbi, Florian Goertz, Diego Redigolo, Robert Ziegler, and Jure Zupan.  Minimal axion model from flavor.  Phys. Rev. D, 95(9):095009, 2017.
  * [69] Sidney R. Coleman and Erick J. Weinberg.  Radiative Corrections as the Origin of Spontaneous Symmetry Breaking.  Phys. Rev. D, 7:1888–1910, 1973.
  * [70] Erick J. Weinberg and Ai-qun Wu.  UNDERSTANDING COMPLEX PERTURBATIVE EFFECTIVE POTENTIALS.  Phys. Rev. D, 36:2474, 1987.
  * [71] James M. Cline and Pierre-Anthony Lemieux.  Electroweak phase transition in two higgs doublet models.  Physical Review D, 55(6):3873–3881, mar 1997.
  * [72] Carroll L. Wainwright.  CosmoTransitions: Computing Cosmological Phase Transition Temperatures and Bubble Profiles with Multiple Fields.  Comput. Phys. Commun., 183:2006–2013, 2012.
  * [73] Rajesh R. Parwani.  Resummation in a hot scalar field theory.  Phys. Rev. D, 45:4695, 1992.  [Erratum: Phys.Rev.D 48, 5965 (1993)].
  * [74] Dario Buttazzo, Giuseppe Degrassi, Pier Paolo Giardino, Gian F. Giudice, Filippo Sala, Alberto Salvio, and Alessandro Strumia.  Investigating the near-criticality of the Higgs boson.  JHEP, 12:089, 2013.
  * [75] Felix Giese and Thomas Konstandin.  Vacuum stability of Froggatt-Nielsen models.  JHEP, 12:091, 2019.
  * [76] Tetsutaro Higaki and Junichiro Kawamura.  A low-scale flavon model with a $Z_{N}$ symmetry.  JHEP, 03:129, 2020.
  * [77] Andrei D. Linde.  Fate of the False Vacuum at Finite Temperature: Theory and Applications.  Phys. Lett. B, 100:37–40, 1981.
  * [78] Archil Kobakhidze, Cyril Lagger, Adrian Manning, and Jason Yue.  Gravitational waves from a supercooled electroweak phase transition and their detection with pulsar timing arrays.  The European Physical Journal C, 77(8), aug 2017.
  * [79] M. Carena, A. Megevand, M. Quirós, and C.E.M. Wagner.  Electroweak baryogenesis and new TeV fermions.  Nuclear Physics B, 716(1-2):319–351, jun 2005.
  * [80] M. Fairbairn and P. Grothaus.  Baryogenesis and dark matter with vector-like fermions.  Journal of High Energy Physics, 2013(10), oct 2013.
  * [81] Hooman Davoudiasl, Ian Lewis, and Eduardo Pontó n.  Electroweak phase transition, higgs diphoton rate, and new heavy fermions.  Physical Review D, 87(9), may 2013.
  * [82] Daniel Egana-Ugrinovic.  The minimal fermionic model of electroweak baryogenesis.  Journal of High Energy Physics, 2017(12), dec 2017.
  * [83] Andrei Angelescu and Peisi Huang.  Multistep strongly first order phase transitions from new fermions at the TeV scale.  Physical Review D, 99(5), mar 2019.
  * [84] Qing-Hong Cao, Katsuya Hashino, Xu-Xiang Li, Zhe Ren, and Jiang-Hao Yu.  Electroweak phase transition triggered by fermion sector.  JHEP, 01:001, 2022.
  * [85] Christophe Grojean, Geraldine Servant, and James D. Wells.  First-order electroweak phase transition in the standard model with a low cutoff.  Phys. Rev. D, 71:036001, 2005.
  * [86] Jose R. Espinosa, Thomas Konstandin, Jose M. No, and Geraldine Servant.  Energy Budget of Cosmological First-order Phase Transitions.  JCAP, 06:028, 2010.
  * [87] Dietrich Bodeker and Guy D. Moore.  Electroweak Bubble Wall Speed Limit.  JCAP, 05:025, 2017.
  * [88] Mark Hindmarsh, Stephan J. Huber, Kari Rummukainen, and David J. Weir.  Numerical simulations of acoustically generated gravitational waves at a first order phase transition.  Phys. Rev. D, 92(12):123009, 2015.
  * [89] Eric Thrane and Joseph D. Romano.  Sensitivity curves for searches for gravitational-wave backgrounds.  Phys. Rev. D, 88(12):124032, 2013.