# Baryogenesis and gravitational waves in the Zee-Babu model

Vo Quoc Phonga,b<EMAIL_ADDRESS>a Department of Theoretical Physics,
University of Science, Ho Chi Minh City 700000, Vietnam
b Vietnam National University, Ho Chi Minh City 700000, Vietnam  Nguyen Chi
Thaoa,b<EMAIL_ADDRESS>a Department of Theoretical Physics,
University of Science, Ho Chi Minh City 700000, Vietnam
b Vietnam National University, Ho Chi Minh City 700000, Vietnam  Hoang Ngoc
Longc,d,e111Corresponding author<EMAIL_ADDRESS>c Subatomic Physics
Group, Science and Technology Advanced Institute, Van Lang University, Ho Chi
Minh City 700000, Vietnam
d Faculty of Technology, Van Lang University University, Ho Chi Minh City
700000, Vietnam
e Institute of Physics, Vietnam Academy of Science and Technology, 10 Dao Tan,
Ba Dinh, Hanoi 100000, Vietnam

###### Abstract

To explain the matter-antimatter asymmetry in the Zee-Babu (ZB) model, the
sphaleron process in the baryogenesis scenario is calculated. It always
satisfies the de-coupling condition and the strength of phase transition ($S$)
is always greater than $1$ in the presence of triggers other than that in the
Standard Model, which are singly ($h^{\pm}$) and doubly ($k^{\pm\pm}$) charged
scalar bosons. Sphaleron energies are in the range of 5-10 TeV, in calculation
with bubble profiles containing free parameters and assuming nuclear bubbles
of $h^{\pm}$ and $k^{\pm\pm}$ are very small. We tested the scaling law of
sphaleron again with an average error of $10\%$. When the temperature is close
to the critical one ($T_{c}$), the density of nuclear bubble is produced very
large and decreases as the temperature decreases. The key parameter is
$\alpha$ which results in the gravitational wave density parameter ($\Omega
h^{2}$) in the range of $10^{-14}$ to $10^{-12}$ when $\beta/H^{*}=22.5$, this
is not enough to detect gravitational waves from electroweak phase transition
(EWPT) according to the LISA data. As the larger strength of phase transition
is, the more $\alpha$ increases (this increase is almost linear with $S$), the
larger the gravitational wave density parameter is. Also in the context of
considering the generation of gravitational waves, in the ZB model we
calculated $\alpha\sim\text{a few}\times 10^{-2}\ll 1$, so rigorously conclude
that the EWPT is not strong even though $S>1$. We also suggest that, for a
model with a lot of extra scalar particles and particles which play a role in
mass generation, the stronger the EWPT process and the larger $\Omega h^{2}$
can be.

###### pacs:

11.15.Ex, 12.60.Fr, 98.80.Cq

Keywords: Spontaneous breaking of gauge symmetries, Extensions of electroweak
Higgs sector, Particle-theory models (Early Universe)

###### Contents

  1. I Introduction
  2. II Review on the Zee-Babu model and the first-order EWPT
  3. III Sphaleron in the Zee-Babu model and proving the sphaleron scaling law
    1. III.1 The zero temperature sphaleron at tree level
    2. III.2 Sphaleron at the non-zero temperature and the scaling law
      1. III.2.1 Sphaleron at the non-zero temperature
      2. III.2.2 The scaling law
  4. IV GW production from EWPT in the Zee-Babu model
    1. IV.1 The parameters $\beta/H^{*}$ and $\alpha$
      1. IV.1.1 The transition rate
      2. IV.1.2 The parameter $\beta/H^{*}$
      3. IV.1.3 The parameter $\alpha$
    2. IV.2 The GW energy density parameter
  5. V Conclusion and discussion

##  I Introduction

The Standard Model (SM) has set up numerous great comes about that concur with
tests and gave us a clear system of how matter interatomic with each other.
The foundation of SM is the concept of symmetries such as Lorentz and gauge
symmetries. The curious thing is that the gauge symmetry unavoidably leads to
the presence of gauge bosons which intercede the intuitive between matter
(fermions). Initially, gauge bosons are massless which are not suitable for
weak interaction acting in short range. Within the Higgs component, which is
fundamentally changes for the situation, can produce masses suddenly,
simultaneously keeps the SM renormalizability.

However, there are a few crucial issues that have not been unraveled by the SM
and one of them is the matter-antimatter asymmetry perplex and we need the
hypothesis to be steady with gravity. We need a dynamical clarification for
this and have to examine electroweak stage prepare that happens within the
early universe. To explain this problem, besides many others, two standard
scenarios are leptogenesis and baryogenesis. In the case of leptogenesis,
according to Ref. 72 , the typical scale is much higher than electroweak
scale. In 1967, Sakharov proposed three conditions that a model must have in
arrange to illuminate the baryogenesis issue sakharov , which are baryon
number violation, C and CP violation (via sphaleron process), and out-of-
equilibrium condition, ie., the first order EWPT and the last condition should
be considered first. In any case, EWPT inside the setting of the SM does not
offer great arrangements since it does not fulfill adequately the two final
conditions of baryogenesis. Hence, models beyond SM got to be taken into
consideration. A common choice is to treat the SM as a complete symmetric
foundation and add only extraneous particles. Of course, we need to come up
with good reasons for adding new particles.

In addition, since the Higgs boson was found at the LHC, Physics Particle
nearly completed its mission that give a more precise understanding of mass.
The EWPT issue will be one of the pressing issues, has been calculated in the
SM as in Refs.mkn ; SME ; michela  and in beyond SM or in numerous other
settings as in Refs.2 ; BSM ; majorana ; thdm ; ESMCO ; elptdm ; phonglongvan
; phonglongvan2 ; SMS ; dssm ; munusm ; lr ; singlet ; mssm1 ; twostep ;
1101.4665 ; jjg ; Ahriche1 ; Ahriche2 ; Ahriche3 ; Fuyuto ; span ; chr ; cde ;
kusenko . We summarize most generally that the triggers for the EWPT process
are either particles outside the SM or manually added quantities 2 ; BSM ;
majorana ; thdm ; ESMCO ; elptdm ; phonglongvan ; phonglongvan2 ; SMS ; dssm ;
munusm ; lr ; singlet ; mssm1 ; twostep ; 1101.4665 ; jjg ; Ahriche1 ;
Ahriche2 ; Ahriche3 ; Fuyuto ; span ; chr ; cde ; kusenko . A commonplace
negative solution is that the EWPT’s strength satisfies the decoupling
condition at a few hundred GeV scale with the mass of Higgs boson must be less
than $125$ GeV mkn ; SME ; michela . So far, the nature of a first-order EWPT
can be non-SM-like particles 2 ; majorana ; thdm ; ESMCO ; elptdm ;
phonglongvan ; phonglongvan2 ; epjc ; zb ; singlet ; mssm1 ; twostep ; chiang3
. Another curiously finding is that EWPT is autonomous of the gauge zb ;
1101.4665 ; Arefe . The damping impact within the self-energy term or daisy
loops have little contributions. When considering daisy loops, the change of
effective potential leads to a diminishment of the strength of EWPT generally
by 2/3 r23 . Subsequently, the contributions do not make an enormous alter to
the strength of EWPT or, in other words, it is not the origin of EWPT.

By ensuring existence C and CP violations, in the context of topological
transitions, the third condition of baryogenesis is defined as
$\Gamma_{sph}\sim\mathcal{O}(T)\textrm{exp}\left(\frac{-E_{sph}}{T}\right)\ll
H_{rad}$ spha-huble ; decoupling , where $\Gamma_{sph}$ and $E_{sph}$ are the
sphaleron rate and energy, respectively, $H_{rad}$ is the Hubble expansion
rate in the radiation-dominated period, and $\mathcal{O}(T)\sim T^{4}$. This
is known as the decoupling condition of sphaleron. In the SM, this condition
is often expressed as $v_{c}/T_{c}>1$ by using an approximation
$E_{sph}(T)\approx[v(T)/v]E_{sph}(T=0)$ 47 ; Ahriche1 ; Ahriche2 ; however, in
beyond SM models this approximation should be again considered with caution.

Sphaleron has also been calculated in various frameworks. We see there are two
main directions: the first is to investigate the phenomena of sphaleron in
models phu1spha ; phu2spha ; phu3spha ; phu4spha ; 10 ; 11 ; 13 ; hagi ; 47 ;
Ahriche1 ; Ahriche2 ; gauge0 , the second is to investigate by simulation and
numerical methods michela ; cosmotran1 . We found that there are two basic
ways to calculate the sphaleron energy plk . The first way, from the sphaleron
energy functional with spherically symmetric ansatzs, deriving the equations
of motion (EOMs), and solving these EOMs and substituting them into the
sphaleron energy functional, we will get desired results. The second way,
assuming the profile functions in the spherically symmetric ansatzs, and
substituting them into the energy functionals, then we will get desired
results. In the first way, we have to use the numerical method twice and have
to check the convergence of the energy functional. In the second way, we only
use the numerical method once, but assuming the profile functions must be very
careful so that the energy functionals converge.

Another aspect, in General Relativity, stars, black holes …are highly
concentrated masses of matter and there are supermassive black hole merger
events, or collisions between two black holes. These collisions are so strong
that they create distortions in space-time, resulting in their energy
densities also changing, which changes the state of space-time and travels out
into space at the speed of light. The propagation of curvature is called a
gravitational wave 1gw . In addition, gravitational waves can be generated by
core collapse of supernovae, the first order EWPT, oscillations of cosmic
strings, inflation. In particular, the strong first order EWPT is one of the
most promising of all these sources 2gw .

There have been many studies on gravitational waves in models phu1gw ; phu2gw
; phu3gw ; phu4gw ; phu5gw ; phu6gwsimulation ; phu7gw ; phu8gw ; 2gw ; 3gw ;
4gw ; 7gw ; 8gw ; 9gw  as well as simulations phu6gwsimulation . Most models
give the values of gravitational wave density parameter smaller than the
observable value (about $10^{-10}$ 7gw ). However, we see that the value of
this parameter may be consistent with the observation if the latent heat
parameter is about $1$.

Among the extended SM models, there are models closest to SM which are The Two
Double Higgs model, SMEFT, the left-right model,… But the ZB model is pretty
simple because the symmetrical structure and particle sectors are like the SM.

In the ZB model, two extra charged scalars $h^{\pm}$ and $k^{\pm\pm}$ are
added to the Higgs content. The kind of new scalars can play an important role
in the early Universe. As shown in Refs.zeebabu ; bzee , they can also be a
reason for tiny mass of neutrinos through two loop or three-loop corrections.

The EWPT in the ZB model was investigated by us in Ref.zb . In the mentioned
paper, the strength of electroweak phase transition ($S$) is greater than one.
We initially concluded that the EWPT in the model under consideration is
strong. However, this is not complete.

To complete the baryogenesis scenario and fully evaluate the EWPT process with
the addition of gravitational waves, are the main purposes of this paper,
which is organized as follows. In section II, we summarize the ZB model and
the EWPT with the one-loop effective potential. In section III, electroweak
sphalerons are calculated at zero and high temperature and we derive again the
sphaleron scaling law. In section IV, $\alpha$ and $\Omega h^{2}$ is
calculated with different strength of phase transition. Finally, we summarize
and make outlooks in section V.

##  II Review on the Zee-Babu model and the first-order EWPT

In the ZB model, by adding two charged scalar bosons $h^{\pm}$ and
$k^{\pm\pm}$ zeebabu , the Lagrangian becomes

$\displaystyle\mathcal{L}$ $\displaystyle=$ $\displaystyle
L_{SM}+f_{ab}\overline{\psi_{aL}^{c}}\psi_{bL}h^{+}+h_{ab}^{{}^{\prime}}\overline{l_{aR}^{c}}l_{bR}k^{++}+V(\phi,h,k)$
(1)
$\displaystyle+(D_{\mu}h^{+})^{\dagger}(D^{\mu}h^{+})+(D_{\mu}k^{++})^{\dagger}(D^{\mu}k^{++})+H.c$

The most important component is the Higgs sector which has four couplings
between $h^{\pm}$ or $k^{\pm\pm}$ and SM-like Higgs zeebabu :

$\displaystyle V(\phi,h,k)$ $\displaystyle=$
$\displaystyle\mu^{2}\phi^{\dagger}\phi+u^{2}_{1}{|{h}|}^{2}+u_{2}^{2}{|{k}|}^{2}+\lambda(\phi^{\dagger}\phi)^{2}+\lambda_{h}{|{h}|}^{4}+\lambda_{k}{|{k}|}^{4}$
(2)
$\displaystyle+\lambda_{hk}{|{h}|}^{2}{|{k}|}^{2}+2p^{2}{|{h}|}^{2}\phi^{\dagger}\phi+2q^{2}{|{k}|}^{2}\phi^{\dagger}\phi+(\mu_{hk}h^{2}k^{++}+H.c)\,,$

where

$\phi=\left(\begin{array}[]{c}\rho^{+}\\\ \rho^{0}\end{array}\right)$ (3)

and $v_{0}$ is a Vacuum Expectation Value (VEV) of $\rho^{0}$,

$\rho^{0}=\frac{1}{\sqrt{2}}\left(v_{0}+\sigma+i\zeta\right)\,.$ (4)

The masses of $h^{\pm}$ and $k^{{\pm\pm}}$ have two components, one comes from
the interaction with the SM-like Higgs boson, the others are free components
($u_{1}^{2}|h|^{2}$ or $u_{2}^{2}|k|^{2}$),

$\displaystyle M^{2}_{h^{\pm}}$ $\displaystyle=$ $\displaystyle
p^{2}v^{2}_{0}+u^{2}_{1}=m^{2}_{h^{\pm}}+u^{2}_{1},m^{2}_{h^{\pm}}=p^{2}v^{2}_{0}\,,$
$\displaystyle M^{2}_{{k^{\pm\pm}}}$ $\displaystyle=$ $\displaystyle
q^{2}v^{2}_{0}+u^{2}_{2}=m^{2}_{{k^{\pm\pm}}}+u^{2}_{2},m^{2}_{{k^{\pm\pm}}}=q^{2}v^{2}_{0}\,.$
(5)

Diagonalizing matrices in the kinetic component of the Higgs potential and
neglect Goldstone bosons, we obtain

$\displaystyle\begin{aligned} m_{H}^{2}(v_{0})&=2\lambda
v^{2}_{0}\,,\end{aligned}\begin{aligned}
m_{Z}^{2}(v_{0})&=\textstyle\frac{1}{4}(g^{2}+g^{\prime
2})v^{2}_{0}=A^{2}v^{2}_{0}\,,\end{aligned}\begin{aligned}
m_{W}^{2}(v_{0})&=\textstyle\frac{1}{4}g^{2}v^{2}_{0}=B^{2}v^{2}_{0}.\end{aligned}$
(6)

The components $u_{1}$ and $u_{2}$ in Eq.(II) do not participate in the EWPT
or the effective potential because they are not components which are derived
from the interaction between Higgs and $h^{\pm}$ or $k^{\pm\pm}$ bosons. This
can be seen clearly in Ref.phonglongvan2 . Therefore, notice that in the next
sections, the masses of $h^{\pm}$ and $k^{\pm\pm}$ appearing hereafter are
just $m_{h^{\pm}}$ and $m_{k^{\pm\pm}}$.

There are usually two ways to calculate the effective potential, which is well
summarized in Ref.plk . We recall the form of the one-loop effective potential
in the Landau gauge zb :

$\displaystyle V_{eff}(v)$ $\displaystyle=$ $\displaystyle
V_{0}(v)+\frac{3}{64\pi^{2}}\left(m_{Z}^{4}(v)ln\frac{m_{Z}^{2}(v)}{Q^{2}}+2m_{W}^{4}(v)ln\frac{m_{W}^{2}(v)}{Q^{2}}-4m_{t}^{4}(v)ln\frac{m_{t}^{2}(v)}{Q^{2}}\right)$
$\displaystyle+$
$\displaystyle\frac{1}{64\pi^{2}}\left(2m_{h^{\pm}}^{4}(v)ln\frac{m_{h^{\pm}}^{2}(v)}{Q^{2}}+2m_{k^{\pm\pm}}^{4}(v){ln}\frac{m_{k^{\pm\pm}}^{2}(v)}{Q^{2}}+m_{H}^{4}(v)ln\frac{m_{H}^{2}(v)}{Q^{2}}\right)$
$\displaystyle+$
$\displaystyle\frac{3T^{4}}{4\pi^{2}}\left\\{F_{-}(\frac{m_{Z}(v)}{T})+F_{-}(\frac{m_{W}(v)}{T})+4F_{+}(\frac{m_{t}(v)}{T})\right\\}$
$\displaystyle+$
$\displaystyle\frac{T^{4}}{4\pi^{2}}\left\\{2{F}_{-}(\frac{m_{h^{\pm}}(v)}{T})+2{F}_{-}(\frac{m_{k^{\pm\pm}}(v)}{T})+{F}_{-}(\frac{{m}_{H}(v)}{T})\right\\}\,,$

where $v$ is the variable that changes with temperature, and at $0K$ then
$v\equiv v_{0}=246$ GeV. $F_{\pm}$ are generalized integrals:

$\displaystyle F_{\pm}\left(\frac{m_{\phi}}{T}\right)$ $\displaystyle=$
$\displaystyle\int_{0}^{\frac{m_{\phi}}{T}}\alpha
J_{\mp}^{1}(\alpha,0)d\alpha,$ $\displaystyle J_{\mp}^{1}(\alpha,0)$
$\displaystyle=$ $\displaystyle
2\int_{\alpha}^{\infty}\frac{(x^{2}-\alpha^{2})^{\frac{1}{2}}}{e^{x}\mp 1}dx.$

After full expansion, with an approximation of $F_{\pm}$, we get:

$V_{eff}(v)=D(T^{2}-T^{2}_{0})v^{2}-ET|v|^{3}+\frac{\lambda_{T}}{4}v^{4}\,.$
(7)

The coefficients in the above potential have the following form:

$\displaystyle D$ $\displaystyle=$
$\displaystyle\frac{1}{24{v_{0}}^{2}}\left[6m_{W}^{2}(v_{0})+3m_{Z_{1}}^{2}(v_{0})+m_{H}^{2}(v_{0})+2m_{h^{\pm}}^{2}(v_{0})+2m_{k^{\pm\pm}}^{2}(v_{0})+6m_{t}^{2}(v_{0})\right],$
$\displaystyle T_{0}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{D}\left\\{\frac{m_{H}^{2}(v_{0})}{4}-\frac{1}{32\pi^{2}v_{0}^{2}}\left(6m_{W}^{4}(v_{0})+3m_{Z_{1}}^{4}(v_{0})+m_{H}^{4}(v_{0})\right.\right.$
$\displaystyle\left.\left.\qquad\qquad\qquad\qquad+2m_{h^{\pm}}^{4}(v_{0})+2m_{k^{\pm\pm}}^{4}(v_{0})-12m_{t}^{4}(v_{0})\right)\right\\},$
$\displaystyle E$ $\displaystyle=$ $\displaystyle\frac{1}{12\pi
v_{0}^{3}}\left(6m_{W}^{3}(v_{0})+3m_{Z_{1}}^{3}(v_{0})+m_{H}^{3}(v_{0})+2m_{h^{\pm}}^{3}(v_{0})+2m_{k^{\pm\pm}}^{3}(v_{0})\right),$
(8) $\displaystyle\lambda_{T}$ $\displaystyle=$
$\displaystyle\frac{m_{H}^{2}(v_{0})}{2v_{0}^{2}}\left\\{1-\frac{1}{8\pi^{2}v_{0}^{2}(m_{H}^{2}(v_{0}))}\left[6m_{W}^{4}(v_{0})\ln\frac{m_{W}^{2}(v_{0})}{a_{b}T^{2}}\right.\right.$
$\displaystyle\qquad\left.\left.+3m_{Z_{1}}^{4}(v_{0})\ln\frac{m_{Z_{1}}^{2}(v_{0})}{a_{b}T^{2}}+m_{H^{0}}^{4}(v_{0})\ln\frac{m_{H^{0}}^{2}(v_{0})}{a_{b}T^{2}}\right.\right.$
$\displaystyle\qquad\left.\left.+2m_{k^{\pm\pm}}^{4}(v_{0})\ln\frac{m_{k^{\pm\pm}}^{2}(v_{0})}{a_{b}T^{2}}+2m_{h^{\pm}}^{4}(v_{0})\ln\frac{m_{h^{\pm}}^{2}(v_{0})}{a_{b}T^{2}}-12m_{t}^{4}(v_{0})\ln\frac{m_{t}^{2}(v_{0})}{a_{F}T^{2}}\right]\right\\}.$

The potential in Eq.(7) has two minima:

$v=0,\quad v\equiv v_{c}=\frac{2ET_{c}}{\lambda_{T_{c}}},$ (9)

where $v_{c}$ is VEV at the critical temperature ($T_{c}$), which is defined:

$T_{c}=\frac{T_{0}}{\sqrt{1-E^{2}/D\lambda_{T_{c}}}}=\frac{T_{0}}{\sqrt{1-ES/2D}}.$
(10)

As the temperature decreases, a second minimum appears. The temperature
$T_{1}$ is the temperature at which the potential begins to appear a second
minimum. $T_{1}$ is calculated as follows

$T_{1}=\frac{T_{0}}{\sqrt{1-9E^{2}/8D\lambda_{T_{c}}}}.$ (11)

The strength of EWPT ($S$) is defined as follows

$S=\frac{v_{c}}{T_{c}}=\frac{2E}{\lambda_{T_{c}}}.$ (12)

As in Ref.zb , $m_{h^{\pm}}$ and beads $m_{k^{\pm\pm}}$ are in the range from
$0$ to $300\,\mathrm{GeV}$, phase transition intensity $1\leq S<2.4$ in the
Landau gauge. However, in that paper, the strength of phase transition can be
up to about 4 with any gauge. Nevertheless, to ensure that the effective
potential does not differ by more than 5%, the mass of the particles must
satisfy the condition $\frac{m}{T}<2.2$ 5percent . Therefore, according to our
full estimation, $S$ cannot increase to 4, is also in the range of 2.4.

We also note that this model has only one mass generation particle (Higgs), so
when $m_{h}$ and $m_{k}$ are too large, $T^{2}_{0}$ in Eq.(8) would be complex
numbers, which makes no sense.

So it can be seen from this model that we want a model with a larger the
strength of phase transition, in addition to having many triggers, the model
must also need more than one SM-like Higgs particle.

##  III Sphaleron in the Zee-Babu model and proving the sphaleron scaling law

In a chiral hypothesis, for example, the electroweak theory there is another
fascinating sort of progress called topological transition, which is a
transition progress between vacua. At the point when such a transition
happens, the baryon number will be violated and it opens up the opportunities
for the electroweak hypothesis to fulfill the main condition of Sakharov. In
this transition, there are two types: Instanton and Sphaleron.

However, with a first order EWPT at high temperature the sphaleron process is
stronger than instanton. The universe is accelerated while this phase
transition occurs thought out sphaleron, so the condition for the departure
from thermal equilibrium (ie., EWPT) is $\Gamma_{sph}\ll H_{rad}$ or the
decoupling condition.

As mentioned in the introduction, there are two ways to calculate the
sphaleron energy. We will do this calculation using the second method
suggested by Manton and Klinkhamer 10 . The ansatz approaches with profiles
functions in Ref.10  to approximate a sphaleron solution without solving the
EOMs.

###  III.1 The zero temperature sphaleron at tree level

As a first step, we have to calculate the sphaleron energy at 0K, which is the
maximum boundary value and as well at the same time to test the convergence of
this functional.

In the ZB model, the sphaleron energy functional consists of three main
components: The first component is gauge bosons, the second and third
components are kinetic and potential of Higgs scalars, respectively,

$E_{sph}=\int
dx^{3}\left[\frac{1}{4}W_{ij}^{a}W_{ij}^{a}+(D_{i}\phi)^{\dagger}(D_{i}\phi)+V(\phi,h^{\pm},k^{\pm\pm})_{tree}\right].$
(13)

We ignore the contributions of the hypercharge gauge field because they are
very small 10 . We also assume the static field approximation (ie.,
$W^{a}_{0}=0$) and sphaleron has a spherically symmetric form 11 :

$\left\\{\begin{array}[]{ll}\phi(r)=\frac{v_{0}}{\sqrt{2}}h(r)in_{a}\sigma\left(\begin{array}[]{cc}0\\\
1\end{array}\right),\\\
W_{i}^{a}(r)=\frac{2}{g}\epsilon^{aij}n_{j}\frac{f(r)}{r},\end{array}\right.$
(14)

where $n_{i}\equiv\frac{x^{i}}{r}$ and $r$ is the radial coordinate in the
spherical coordinates. Here, $h(r)$ and $f(r)$ have boundary conditions

$\left\\{\begin{array}[]{ll}h(r\to 0)=f(r\to 0)=0,\\\
h(r\to\infty)=f(r\to\infty)=1.\end{array}\right.$ (15)

There are two types of profile ansatzs: the ansatz with scale-free parameters
firstly introduced in Ref. 10 , and the smooth ansatz as in Ref. 11 ; gauge0 .
In addition, with composite Higgs, models can have other profiles phu2spha .
According to the calculation in Ref.plk , the ansatz with scale-free
parameters is more accurate. So we use the ansatz with two scale-free
parameters $a$ and $b$:

$f(\xi)=\left\\{\begin{array}[]{c}f_{1}(\xi)=\frac{\xi^{2}}{2a^{2}},\xi\leq
a\\\ f_{2}(\xi)=1-\frac{a^{2}}{\xi^{2}},\xi\geq a\par\end{array}\right.$ (16)
$h(\xi)=\left\\{\begin{array}[]{c}h_{1}(\xi)=\frac{4\xi}{5b},\xi\leq b\\\
h_{2}(\xi)=1-\frac{b^{4}}{5\xi^{4}},\xi\geq b\end{array}\right.$ (17)

where

$\xi\equiv gv_{0}r\rightarrow d\xi=gv_{0}dr.$

The potential at tree level can be rewritten as follows:

$\displaystyle V(\phi,h^{\pm},k^{\pm\pm},T=0)_{tree}$ $\displaystyle=$
$\displaystyle\mu^{2}\phi^{\dagger}\phi+\lambda(\phi^{\dagger}\phi)^{2}$ (18)
$\displaystyle+$ $\displaystyle
2p^{2}{|{h}|}^{2}\phi^{\dagger}\phi+2q^{2}{|{k}|}^{2}\phi^{\dagger}\phi+\mathcal{V}\,$

where

$\displaystyle\mathcal{V}$ $\displaystyle=$ $\displaystyle
u^{2}_{1}{|{h}|}^{2}+u_{2}^{2}{|{k}|}^{2}+\lambda_{h}{|{h}|}^{4}+\lambda_{k}{|{k}|}^{4}+\lambda_{hk}{|{h}|}^{2}{|{k}|}^{2}$
(19) $\displaystyle+$ $\displaystyle(\mu_{hk}h^{2}k^{++}+H.c).$

We next accept an approximation:

$\displaystyle h^{\pm}(r)$ $\displaystyle=$ $\displaystyle c_{1},$ (20)
$\displaystyle k^{\pm\pm}(r)$ $\displaystyle=$ $\displaystyle c_{2},$ (21)

where $c_{1}$ and $c_{2}$ are very small. So $\mathcal{V}$ in Eq. (19) becomes
constant which can be set to $\Lambda$.

We were able to make this approximation because the bubble nucleation is the
process of mass generation of particles through interaction with Higgs field.
$\mathcal{V}$ contains terms that do not interact to Higgs boson. This is
approximate, but may be required because $\mathcal{V}$ is entered manually and
they must satisfy a few conditions, so as not to disrupt the natural
properties of the model. So if the bubbles of them are too large, it will make
$\mathcal{V}$ to have a very large contribution to the sphaleron energy,
possibly a dominant component, which makes no sense. Also, the coefficients in
$\mathcal{V}$ are undefined. The values of these coupling constants need to be
determined through decay channels. Therefore, in order to have the keys to
handle other problems and not lose the naturalness of the model, our
approximation is feasible.

After putting the ansatzs into Eq. (18), we will get

$\displaystyle V(\phi,h^{\pm},k^{\pm\pm},T=0)_{tree}$ $\displaystyle=$
$\displaystyle-\frac{\mu^{\prime 2}v_{0}^{2}}{2}h^{2}+\frac{\lambda
v_{0}^{4}}{4}h^{4}+\Lambda$ (22)

where $\mu^{\prime
2}=\mu^{2}+4p^{2}c^{2}_{1}+4q^{2}c^{2}_{2}=constant\sim\mu^{2}$. We can
minimize the above potential and eliminate $\Lambda$, then obtain

$\displaystyle V(\phi,h^{\pm},k^{\pm\pm},T=0)_{tree}=\frac{\lambda
v_{0}^{4}}{4}(h^{2}-1)^{2}.$ (23)

After making a few simple transformations, the zero temperature sphaleron
energy at tree-level is

$E_{sph,tree}^{0K}=\frac{4\pi
v_{0}}{g}\int_{0}^{\infty}d\xi\left[4\left(\frac{df}{\xi}\right)^{2}+\frac{8f^{2}(1-f)^{2}}{\xi^{2}}+\frac{\xi^{2}}{2}\left(\frac{dh}{\xi}\right)^{2}+h^{2}(1-f)^{2}+\frac{\xi^{2}}{g^{2}}\frac{\lambda}{4}(h^{2}-1)^{2}\right]$
(24)

Expanding Eq.(24) by using Eqs. (16) and (17), we get

$\displaystyle E_{sph,tree}^{0K}$ $\displaystyle=$ $\displaystyle\frac{4\pi
v_{0}}{g}\int_{0}^{a}d\xi\left[4\left(\frac{df_{1}(\xi)}{\xi}\right)^{2}+\frac{8f^{2}_{1}(\xi)(1-f_{1}(\xi))^{2}}{\xi^{2}}\right]$
(25) $\displaystyle+$ $\displaystyle\frac{4\pi
v_{0}}{g}\int_{0}^{a}d\xi\left[\frac{\xi^{2}}{2}\left(\frac{df_{1}(\xi)}{\xi}\right)^{2}+\frac{8h^{2}_{1}(\xi)(1-f_{1}(\xi))^{2}}{\xi^{2}}+\frac{\xi^{2}\lambda}{4g^{2}}\left(h_{1}^{2}(\xi)-1\right)^{2}\right]$
$\displaystyle+$ $\displaystyle\frac{4\pi
v_{0}}{g}\int_{a}^{b}d\xi\left[4\left(\frac{df_{2}(\xi)}{\xi}\right)^{2}+\frac{8f^{2}_{2}(\xi)(1-f_{2}(\xi))^{2}}{\xi^{2}}\right]$
$\displaystyle+$ $\displaystyle\frac{4\pi
v_{0}}{g}\int_{a}^{b}d\xi\left[\frac{\xi^{2}}{2}\left(\frac{dh_{1}(\xi)}{\xi}\right)^{2}+\frac{8h^{2}_{1}(\xi)(1-f_{2}(\xi))^{2}}{\xi^{2}}+\frac{\xi^{2}\lambda}{4g^{2}}\left(h_{1}^{2}(\xi)-1\right)^{2}\right]$
$\displaystyle+$ $\displaystyle\frac{4\pi
v_{0}}{g}\int_{b}^{\infty}d\xi\left[4\left(\frac{dh_{2}(\xi)}{\xi}\right)^{2}+\frac{8h^{2}_{2}(\xi)(1-h_{2}^{2}(\xi))}{\xi^{2}}\right]$
$\displaystyle+$ $\displaystyle\frac{4\pi
v_{0}}{g}\int_{b}^{\infty}d\xi\left[\frac{\xi^{2}}{2}\left(\frac{dh_{2}(\xi)}{\xi}\right)^{2}+\frac{8h^{2}_{2}(\xi)(1-f_{2}(\xi))^{2}}{\xi^{2}}+\frac{\xi^{2}\lambda}{4g^{2}}\left(h_{2}^{2}(\xi)-1\right)^{2}\right],$

where $g^{2}=\frac{G_{F}8m^{2}_{W}}{\sqrt{2}}$; $G_{F}=1.166\times
10^{-}{15}\,\textrm{GeV}^{-2}$; $m_{W}=80.39\,\textrm{GeV}$; $\lambda\simeq
0.12$.

To minimize this energy according to the parameters $a$ and $b$, we get the
values of sphaleron energy given in Table 1.

Table 1: The zero temperature sphaleron energy at tree level $a$ | $b$ |  $v_{0}$[GeV] | $E_{sph}^{T=0K}[\textrm{GeV}]$ | $E_{sph}^{T=0K}/v_{0}$
---|---|---|---|---
$2.4$ | $1.28$ | 246 | $10245.7$ | $41.649$

Note that Eq.(25) has a symmetry between $a$ and $b$. If we assume $a>b$ or
$b>a$, the above integral gives the same minimum value. The maximum sphaleron
energy which is one at 0K, can be estimated at around $10.2457$ TeV.

###  III.2 Sphaleron at the non-zero temperature and the scaling law

####  III.2.1 Sphaleron at the non-zero temperature

We just need to replace $V(\phi,h^{\pm},k^{\pm\pm})$ by $V(h,T)_{eff}$ in the
sphaleron energy at 0K, then we will get sphaleron energy at the non-zero
temperature

$E_{sph}^{T}=\frac{4\pi
v_{0}}{g}\int_{0}^{\infty}d\xi\left[4\left(\frac{df}{\xi}\right)^{2}+\frac{8f^{2}(1-f)^{2}}{\xi^{2}}+\frac{\xi^{2}}{2}\left(\frac{dh}{\xi}\right)^{2}+h^{2}(1-f)^{2}+\frac{\xi^{2}}{g^{2}v^{4}_{0}}V_{eff}(h,T)\right],$
(26)

where

$V_{eff}(h,T)=D(T^{2}-T^{2}_{0})v^{2}-ET|v|^{3}+\frac{\lambda_{T}}{4}v^{4}.$
(27)

Notice that

$\phi^{\dagger}\phi=\frac{v_{0}^{2}h^{2}(r)}{2}\,.$ (28)

On the other hand, ignoring the contribution of the Goldstone bosons, the
Higgs doublet becomes

$\phi=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}0\\\
v\end{array}\right)\to\phi^{\dagger}\phi=\frac{v^{2}}{2}.$ (29)

Form (28) and (29), we obtain

$v=v_{0}\times h(r)$ (30)

Inserting (30) into (27) we then get

$\displaystyle
V_{eff}^{T}=D(T^{2}-T^{2}_{0})v_{0}^{2}h^{2}-ET|h|^{3}v_{0}^{3}+\frac{\lambda_{T}}{4}v_{0}^{4}h^{4}.$
(31)

We must rescale the variable $\xi$ and the radial function $h(\xi)$ as
follows:

$\tilde{\xi}=\frac{\xi}{s};\hskip 28.45274pt\tilde{h}=sh;\hskip
28.45274pt\text{with }s\equiv\frac{v(T)}{v_{0}}.$

then Eq. (31) becomes

$\displaystyle V_{eff}(h,T)$ $\displaystyle=$
$\displaystyle\frac{\lambda_{T}}{4}v^{4}-ETv^{3}+D(T^{2}-T_{0}^{2})v^{2}$ (32)
$\displaystyle=$
$\displaystyle\frac{\lambda_{T}v_{0}^{4}}{4}s^{4}h^{4}-ETv_{0}^{3}s^{3}h^{3}+D(T^{2}-T_{0}^{2})v_{0}^{2}s^{2}h^{2}$
$\displaystyle=$
$\displaystyle\frac{\lambda_{T}v_{0}^{4}}{4}s^{4}(h^{4}-1)-ETv_{0}^{3}s^{3}(h^{3}-1)+D(T^{2}-T_{0}^{2})v_{0}^{2}s^{2}(h^{2}-1)$
$\displaystyle+$ $\displaystyle\mathcal{C}(T).$

In Eq.(32), $\mathcal{C}(T)$ will be removed by the normalized conditions of
the effective potential, so that $E_{sph}^{T}$ does not diverge. The sphaleron
energy will be rewritten as

$\displaystyle E_{sph}^{T}$ $\displaystyle=$ $\displaystyle\frac{4\pi
v_{0}}{gs}\int_{0}^{\infty}d\xi\Bigg{\\{}s^{2}\Bigg{[}4\left(\frac{df}{d\xi}\right)^{2}+\frac{8}{\xi^{2}}f^{2}(1-f)^{2}+\frac{\xi^{2}}{2}\left(\frac{dh}{d\xi}\right)^{2}+h^{2}(1-f)^{2}\Bigg{]}$
(33) $\displaystyle+$
$\displaystyle\frac{\xi^{2}}{g^{2}s^{2}}\Bigg{[}\frac{\lambda_{T}}{4}s^{4}(h^{4}-1)-\frac{ET}{v_{0}}s^{3}(h^{3}-1)+\frac{D(T^{2}-T_{0}^{2})}{v_{0}^{2}}s^{2}(h^{2}-1)\Bigg{\\}}.$

In Eq. (33) we can separate out the sphaleron energy of the gauge fields in
terms of the variable $\xi$ as follows plk

$\displaystyle E^{T}_{gauge}$
$\displaystyle=\int\frac{1}{4}W_{ij}^{a}W_{ij}^{a}dx^{3}=\frac{4\pi
v}{g}\int_{0}^{\infty}d\xi\left[4\left(\frac{df}{d\xi}\right)^{2}+\frac{8f^{2}(1-f)^{2}}{\xi^{2}}\right].$
(34)

With different $S$ values, we get the specific mass of $h^{\pm}$ and
$k^{\pm\pm}$. Then substitute them in Eq.(33), and minimize the energy
functional to find the values $a,b$. Finally, substituting $a,b$ which we just
found with the mass of the particles in Eq.(33), then we get the sphaleron
energy at any temperature in the range below $T_{c}$.

To ensure that the baryon number remains unchanged during the expansion of the
universe, we need the baryon erosion avoidance condition, which is also called
the de-coupling condition or the third Sakharov condition. During the phase
transition (that is, when the temperature drops from $T_{c}$ to $T_{0}$), the
rate of violation of the baryon number must be less than the Hubble rate.
Moreover, EWPT occurs after inflation and belongs to the period dominated by
radiation, so the Hubble rate is
$H^{2}_{rad}=\frac{4\pi^{3}g_{*}}{45m_{pl}^{2}}T^{4}$. Thus $\Gamma_{sph}\ll
H_{rad}$ spha-huble ; decoupling ; 11  as the result:

$\frac{E_{sph}^{T}}{T}-7ln\left(\frac{v(T)}{T}\right)+ln\left(\frac{T}{100\,\textrm{GeV}}\right)>(35.9-42.8).$
(35)

For detailed illustration, we calculate the sphaleron energy with $S=1.2$, the
results obtained are as shown in Table 2. These results satisfy the decoupling
and supercooling conditions Fuyuto .

Table 2: The sphaleron energy at $S=1.2$, $m_{k^{\pm\pm}}=256.169\,\textrm{GeV}$, $m_{h^{\pm}}=218.500\,\textrm{GeV}$. $T[\textrm{GeV}]$ | $v[\textrm{GeV}]$ | $a$ | $b$ | $E^{T}_{gauge}$ | $E_{sph}^{T}[\textrm{GeV}]$ | $E^{T}_{gauge}/E_{sph}^{T}$ | De-coupling: Eq.(35)
---|---|---|---|---|---|---|---
$T_{c}=118.104$ | 141.724 | 2.591 | 2.855 | 2792.88 | 5245.45 | 53.24% | $43.304>42.8$
$T_{n1}=117.000$ | 150.056 | 2.540 | 2.715 | 3016.16 | 5612.83 | 53.73% | $46.38>42.8$
$T_{n2}=115.000$ | 161.346 | 2.490 | 2.579 | 3307.95 | 6103.09 | 54.20% | $50.83>42.8$
$T_{n3}=113.000$ | 170.055 | 2.461 | 2.499 | 3528.46 | 6478.11 | 54.46% | $54.58>42.8$
$T_{n4}=111.000$ | 177.238 | 2.440 | 2.444 | 3708.61 | 6786.33 | 54.64% | $57.96>42.8$
$T_{0}=105.794$ | 191.66 | 2.420 | 2.425 | 4042.93 | 7407.11 | 54.58% | $65.91>42.8$

We see that $E_{sph}^{T}$ at any temperatures satisfy this condition, see in
the eighth column of the Table 2, the gauge component contributed the largest,
accounting for about 55%. Although the results are only examined for one value
of $S$, the laws of sphaleron energy hold for other $S$.

The profile functions of $h(\xi)$ and $f(\xi)$ are plotted in Fig.1. These
functions rapidly advance to $1$ as $\xi$ increases. As $\xi$ increases, the
faster the functions approach $1$, the thinner the wall of bubbles.

(a) The $h(\xi)$ profile with $S=1.2$. The red line: $T=T_{c}$, the blue and
black line: $T=T_{0}$ and $T_{n2}$.

(b) The $f(\xi)$ profile with $S=1.2$. The red line: $T=T_{c}$, the black
line: $T=T_{0}$.

Figure 1: Nuclear bubbles

According to the results of Table 2, the values of $a,b$ decrease as the
temperature decreases. When the temperature is near $T_{c}$, $a$ and $b$
decrease very quickly (see Fig. 2). $a,b$ are two free parameters, but they
are specific to the thickness of the wall of bubbles. This leads to an initial
prediction that the radius of bubbles will increase when the temperature goes
down because the thickness of wall of bubbles decrease. We will verify this
again in the next section after calculating $\beta/H^{*}$.

Figure 2: The values of a and b depend on $T$

####  III.2.2 The scaling law

The sphaleron scaling law 47 ; Ahriche1 ; Ahriche2 : The sphaleron energy at a
temperature below the critical temperature is linearly proportional to VEV at
that temperature

$E_{scaling}^{T}\approx\frac{v(T)}{v_{0}}E_{sph}^{T=0K}.$ (36)

This law is only an approximate expansion. But it has a necessary meaning in
order to double-check the calculation results.

Using data in Table 2, we have plotted $E_{sph}^{T}$ (in the sixth column) as
a function of $v(T)$ (the second column) in Fig. 3. From Fig. 3, $E^{T}_{sph}$
increases linearly with $v(T)$ and the slope of this line has the typical
range $37-38$. But according to Eq.(36), this slope should be equal to
$\frac{E_{sph}^{0K})}{v_{0}}$, so that
$\tan(\theta)=\frac{E_{sph}^{0K}}{v_{0}}$. From the Table 1, we get
$\tan(\theta)=\frac{E_{sph}^{0K}}{v_{0}}=\frac{10245.7}{246}=41.649$. We
therefore find that the above law is on average 10% more than it actually is.

Figure 3: The sphaleron scaling law

Hence, we revise the sphaleron scaling law as follows:

$E_{scaling}^{T}\approx v(T).\frac{0.9E_{sph}^{0K}}{v_{0}}.$ (37)

This can be clearly seen when comparing Eq.(25) with Eq.(26); and
$V(\phi,h^{\pm},k^{\pm\pm})$ is very different from $V_{eff}(h,T)$. However,
this correction is not large because the contribution of the effective
potential component is small, which is also mentioned in Ref. Ahriche1 . In
addition this correction is also consistent with the results given in Ref.plk
.

##  IV GW production from EWPT in the Zee-Babu model

Electroweak phase transition with nuclear bubbles that has spherical symmetry.
As the temperature drops from the critical temperature, the collisions of
bubbles destroy the spherical symmetry. This produces gravitational radiation.
The high velocity and energy density of the bubbles will provide the necessary
conditions for the generation of gravitational radiation. There are two
sources of gravitational waves: the direct collision of the bubbles and the
chaos in the plasma state due to the motion of the bubbles 2gw ; 7gw .

###  IV.1 The parameters $\beta/H^{*}$ and $\alpha$

Two important parameters for calculating gravitational waves that have been
derived from many past studies phu1gw ; phu2gw ; phu3gw ; phu4gw ; phu5gw ;
phu7gw ; phu8gw ; 1gw ; 2gw ; 3gw ; 4gw ; 7gw ; 8gw ; 9gw  are latent heat
($\alpha$) and $\beta/H^{*}$.

The first parameter, $\beta$ is often called the inverse of the time of the
phase transition. This quantity determines the size of the bubble at the time
of collision, so it is characteristic of the frequency at which the GW signal
is largest 1gw .

####  IV.1.1 The transition rate

The sphaleron rate below the critical temperature is usually defined as:

$\Gamma(t)=A(t)e^{-S(t)}.$ (38)

We deduce the derivative of this rate with respect to time:

$\frac{d\Gamma(t)}{dt}=\frac{dA(t)}{dt}e^{-S(t)}-A(t)\frac{dS(t)}{dt}e^{-S(t)}.$
(39)

in which, $A(t)\sim T^{4}$ so
$\frac{dA(t)}{dt}\approx\frac{d(T^{4})}{dt}=\frac{d(T^{4})}{dT}\frac{dT}{dt}=4T^{3}\frac{dT}{dt}$
and $S(t)$ is the Euclidean action of a critical bubble.

The change of temperature over time, calculated through the Hubble coefficient
$H$ gives $\frac{dT}{dt}=-T\times H$. So (38) becomes:

$\displaystyle\frac{d\Gamma(t)}{dt}=-4T^{4}He^{-S(t)}+T^{5}H\frac{d\left(S(t)\right)}{dT},$
(40)

####  IV.1.2 The parameter $\beta/H^{*}$

The parameter $\beta$ is defined as follows:

$\displaystyle\beta$ $\displaystyle\approx$
$\displaystyle\frac{1}{\Gamma}\frac{d\Gamma(t)}{dt}=\left[-4H+TH\frac{d\left(S(t)\right)}{dT}\right],$
(41)

with $H\approx 3.36\times 10^{-14}$. So we can ignore the term $-4H$ in
Eq.(41) so it can be rewritten as follows

$\displaystyle\frac{\beta}{H^{*}}$ $\displaystyle\approx$
$\displaystyle\left[T\frac{d\left(S(t)\right)}{dT}\right]_{T=T_{C}}.$ (42)

We can approximate $S(t)\approx S_{3}(T)/T$. Here $S_{3}(T)$ is the action for
the $O(3)$ symmetric bounce action showing transitions between vacua of the
Higgs potential, excluding gauge fields. Therefore $S_{3}(T)$ can be
approximated by sphaleron energy without the gauge component. From the values
in Table 2, it follows that the gauge contribution accounts for 55% of the
sphaleron energy. Therefore, a possible approximation is

$\frac{\beta}{H^{*}}\approx\left[T\frac{d\left(S(t)\right)}{dT}\right]_{T=T_{C}}\approx\left[T\frac{d\left(S_{3}/T\right)}{dT}\right]_{T=T_{C}}\approx\left[T\frac{d\left(0.45\frac{E^{T}_{sph}}{T}\right)}{dT}\right]_{T=T_{C}}.$
(43)

Here, $\frac{\beta}{H^{*}}$ is dimensionless and depends mainly on the
sphaleron energy and the shape of the potential at the time of nucleation.
$H^{*}$ is defined as follows 8gw

$\displaystyle H^{*}$ $\displaystyle=$ $\displaystyle 4.5\times
10^{-22}\left[\frac{T_{c}}{MeV}\right]^{2}\,\textrm{
MeV}\left[\frac{g_{*}}{10.75}\right]^{2};\,g^{*}=106.75.$ (44)

In addition we have three quantities that are related to the nuclear bubble
defined via the $\beta$ parameter,

$\displaystyle\begin{cases}\beta^{-1}\Gamma(t):\text{the number of nuclear
bubbles,}\\\ R_{b}\sim v_{b}\beta^{-1}:\text{the radius of nuclear bubble,}\\\
\beta^{-1}:\text{the typical time-scale of transition process.}\end{cases}$
(45)

in which, $v_{b}$ is the velocity of the wall of bubble.

The parameter $\beta$ is difficult to determine, because we need to know the
change of the sphaleron energy functional around the critical temperature. But
Eq.(26) only tells us the energy at $T\leq T_{c}$. And the value of the
sphaleron energy when $T>T_{c}$ is often difficult to determine, because the
sphaleron rate in the thermal equilibrium is very small due to the sphaleron
processes can be wash-outed by the inverse process.

From Eq. (43) it is reasonable the following approximation:

$\beta/H^{*}\approx\bigg{[}\frac{0.45E^{T}_{sph}}{T}\bigg{]}_{T=T_{c}}.$ (46)

####  IV.1.3 The parameter $\alpha$

The second important parameter characterizing the spectrum of gravitational
waves is $\alpha$ that is the ratio of latent heat ($\epsilon$) released in
the radiation energy density ($\rho^{*}_{rad}$). $\alpha$ is determined at the
nucleation temperature as follows 2gw :

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\frac{\epsilon}{\rho_{rad}^{*}},$ $\displaystyle\epsilon$
$\displaystyle=$
$\displaystyle\left(-V(v(T),T)+T\frac{d}{dT}V(v(T),T)\right)_{T=T_{C}},$
$\displaystyle\rho_{rad}^{*}$ $\displaystyle=$ $\displaystyle
g^{*}\pi^{2}\frac{T^{4}_{c}}{30}=106.75\pi^{2}\frac{T^{4}_{c}}{30}.$ (47)

This parameter depends only on the configuration of the effective potential.
Thus in all models the parameter which is easily calculated by usual way when
the effective potential is known, determines the strength of the phase
transition.

###  IV.2 The GW energy density parameter

With the effective potentials in the above sections, we can compute the
quantities $\alpha$ and $\beta/H^{*}$ and from that we can predict for the GW
spectrum which is characterized by the density parameter.

There are three processes which are involved in the production of GWs at a
first-order PT: Collisions of bubble walls, sound waves in the plasma and
Magnetohydrodynamic (MHD) turbulence 71 . These three processes generically
coexist, and the corresponding contributions to the stochastic GW background
should linearly combine, at least approximately, so that 71

$h^{2}\Omega_{GW}=h^{2}\Omega_{Coll}+h^{2}\Omega_{sw}+h^{2}\Omega_{tur}.$ (48)

We have three different scenarios depending on the bubble wall velocity 71 .
The bubble wall velocity is defined via $\alpha$:

$v_{b}=\frac{\frac{1}{\sqrt{3}}+\sqrt{\alpha^{2}+2\frac{\alpha}{3}}}{1+\alpha}.$
(49)

We estimate that $v_{b}\sim 1$. This will be clearly seen on Table 3.
Therefore, the generation of gravitational waves has only two components: the
sound wave and turbulence 71 ; 73 .

The gravitational wave energy density parameters of sound waves is defined
through $\alpha$ and $\beta/H^{*}$ as follows phu5gw ; 71 :

$h^{2}\Omega_{sw}(f)\simeq 2.65\times
10^{-6}v_{b}k^{2}_{w}\left[\frac{H_{*}}{\beta}\right]\left[\frac{\alpha}{1+\alpha}\right]^{2}\left[\frac{100}{g^{*}}\right]^{\frac{1}{3}}S_{sw}(f),$
(50)

where

$S_{sw}(f)=(f/f_{sw})^{1/3}\left(\frac{7}{4+3(f/f_{peak.sw})^{2}}\right)^{7/2},$
(51)

and the peak frequency is well fitted by

$f_{peak.sw}=1.9\times
10^{-2}mHz\times\frac{1}{v_{b}}\frac{\beta}{H^{*}}\frac{T_{*}}{100}\left(\frac{g_{*}}{100}\right)^{1/6}.$
(52)

GWs from turbulence in the cosmic fluid is given by

$h^{2}\Omega_{tur}(f)\simeq 3.35\times
10^{-4}\left(\frac{H^{*}}{\beta}\right)v_{b}\left(\frac{k_{t}\alpha}{1+\alpha}\right)^{\frac{3}{2}}\left(\frac{100}{g^{*}}\right)^{1/3}S_{tur}(f),$
(53)

in which

$S_{tur}(f)=\frac{(f/f_{peak.tur})^{3}}{\left(1+(f/f_{peak.tur})\right)^{11/3}\left(1+8\pi
f/h^{*}\right)},$ (54)

and Eq. (54) reads

$\displaystyle h_{*}$ $\displaystyle=$ $\displaystyle 16.5\times
10^{-3}mHz\times\frac{T_{*}}{100}\left(\frac{g_{*}}{100}\right)^{1/6}\,,$
$\displaystyle f_{peak.tur}$ $\displaystyle=$ $\displaystyle 2.7\times
10^{-2}mHz\times\frac{1}{v_{b}}\frac{\beta}{H^{*}}\frac{T_{*}}{100}\left(\frac{g_{*}}{100}\right)^{1/6}.$
(55)

One can estimate that 71

$\displaystyle
k_{w}=\alpha\left[0.73+0.083\sqrt{\alpha}+\alpha\right]\quad\text{when}\quad
v_{b}\sim 1;$ (56) $\displaystyle k_{t}=(0.05-0.1)\times k_{w}.$ (57)

where $k_{w}$ is part of the latent heat, converted into kinetic energy of the
bubble wall.

It is important to note that the current formulas for the calculation of the
gravitational wave spectrum have some differences, possibly slightly from
Eqs.(50) and (53), but they are all the same power. One of the pretty good
corrections available today as in Ref.grahamgw1 . However, in order for us to
only evaluate the power of the gravitational wave spectrum, we only use the
formulas in Eqs.(50) and (53).

First, we have to evaluate $\alpha$ and $\beta/H^{*}$: if $\alpha$ increases
and $\beta/H^{*}$ decreases, the density parameter increases. This is shown as
Fig.4.

(a) $\Omega h^{2}$ with different $\beta/H^{*}$ values, $f=f_{peak}$

(b) $\Omega h^{2}$ with different $\alpha$ values, $f=f_{peak}$

Figure 4: The graph of $\Omega h^{2}$ varies with $\alpha$ and $\beta/H^{*}$

So, according to Eqs.(43) and (47) if we want the larger density parameter to
be, the larger $\alpha$ is, i.e., the larger change of the potential over
temperature should be, and the smaller the change of the sphaleron energy over
temperature is at the critical temperature.

By approximating Eq.(46) and the data in Table 2, we get $\beta/H^{*}\approx
22.5$. We then calculate $\alpha$ and the values of $\Omega h^{2}$ with a few
$S$, as shown in the Table 3.

Table 3: Gravitational wave corresponding to the strength of phase transition with $\beta/H^{*}=22.5$ and $f=f_{peak}$. $T_{c}[\textrm{GeV}]$ | $v(T_{c})[\textrm{GeV}]$ | $\alpha$ | $[\Omega h^{2}]_{sw}$ | $[\Omega h^{2}]_{tur}$ | $\Omega h^{2}]_{ewpt}$ | $v_{b}$
---|---|---|---|---|---|---
$S=1$, $m_{k^{\pm\pm}}=219.496\,\textrm{GeV}$,
$m_{h^{\pm}}=209.293\,\textrm{GeV}$
$123.081$ | $123.08$ | $0.0148$ | $6.56499\times 10^{-14}$ | $2.93648\times 10^{-14}$ | $3.59298\times 10^{-14}$ | $0.668$
$S=1.2$, $m_{k^{\pm\pm}}=256.169\,\textrm{GeV}$,
$m_{h^{\pm}}=218.500\,\textrm{GeV}$
$118.099$ | $141.724$ | $0.022$ | $3.13011\times 10^{-14}$ | $9.45874\times 10^{-14}$ | $1.25888\times 10^{-13}$ | $0.685$
$S=1.4$, $m_{k^{\pm\pm}}=294.418\,\textrm{GeV}$,
$m_{h^{\pm}}=217.491\,\textrm{GeV}$
$113.606$ | $159.049$ | $0.0304$ | $1.08423\times 10^{-13}$ | $2.39771\times 10^{-13}$ | $3.48193\times 10^{-13}$ | $0.701$
$S=1.6$,$m_{k^{\pm\pm}}=288.500\,\textrm{GeV}$,
$m_{h^{\pm}}=278.355\,\textrm{GeV}$
$108.986$ | $174.379$ | $0.0403$ | $3.26923\times 10^{-13}$ | $5.47717\times 10^{-13}$ | $8.7464\times 10^{-13}$ | $0.717$
$S=1.8$, $m_{k^{\pm\pm}}=318.163\,\textrm{GeV}$,
$m_{h^{\pm}}=287.461\,\textrm{GeV}$
$104.934$ | $188.881$ | $0.0505$ | $7.81185\times 10^{-13}$ | $1.05107\times 10^{-12}$ | $1.83226\times 10^{-12}$ | $0.730$

According to the Table 3, as $S$ increases, $v_{b}$ increases. So the radius
of bubbles expressed in Eq. (45) also increases. We next look at the first and
last columns of the Table 3, as the temperature decreases, $v_{b}$ increases
again, which tells us that as the temperature decreases the radius of bubbles
increases, or the velocity of the walls of bubbles increases, leads to the
thickness of the walls can be reduced. This was also commented earlier, when
observing the values of the two free parameters $a$ and $b$.

As the temperature decreases, from the Tables 2 and 3 with Eq.(45), the
sphaleron energy increases so the sphaleron rate decreases, consequently the
number of bubbles also decreases.

According to the figures 5, we see that $\Omega h^{2}$ increases when $S$ is
greater than $1$ but not yet in the range $10^{-12}-10^{-10}$ and the value of
$\alpha$ increases linearly with $S$. So in the model under consideration, the
phase transition strength should be around $3.6$, $\alpha$ would have a value
in the range of $1$. The last, we conclude that, with the strength of phase
transition of about $3.6$, this model shows that the gravitational wave energy
density is within the observable region of the LISA. But the maximum strength
in the model under consideration is only about $2.4$, so that the error of the
effective potential does not exceed 5% 5percent .

(a) The change of $\alpha$ with the strength of phase transition.

(b) The change of $\Omega h^{2}$ with the strength of phase transition.

Figure 5: Graphs of $\alpha$ and $\Omega h^{2}$, $f=f_{peak}$.

##  V Conclusion and discussion

The sphaleron energy has been calculated by profile functions with two free
parameters. This is a scenario with good results, the numerical solution is
not complicated, but the convergence of the energy functional is fast.
However, we need to carefully deal with the free constants in the energy
functional before solving.

We have noticed that the sphaleron energy needs to satisfy the following,
regardless of the way: the first is the de-coupling condition; the second is
the scaling law; the maximum value of sphaleron energy is at zero temperature;
and the last, we can only compute the sphaleron energy in detail at
temperature $T\leq T_{c}$ and can only estimate the rate of sphaleron at
temperatures above $T_{c}$.

The sphaleron rate tells us about the baryon-number violation, the first
condition of baryogenesis. In the ZB model has the nonzero phases of the CKM
matrix and a neutrino mixing 74 , they tell us about CP violation. The last
condition for thermal imbalance is demonstrated by a first order EWPT in this
model. Thus, the ZB model allows a baryogenesis scenario that explains the
matter-antimatter asymmetry.

As estimated for $\Omega h^{2}$ to be around $10^{-12}-10^{-10}$ (within the
LISA observation limit), $\alpha$ must be greater than $1$, $S$ must be
greater than $3$. However, $S$ cannot be very large (should only be in the
vicinity of $2$), otherwise the effective potential expansions are no longer
accurate. An error can also come from the slow convergence of perturbation
theory in the Quantum Field Theory. This can be greatly improved after the re-
summation process as shown in Refs. grahamgw2 .

Our calculations still satisfy the condition $m/T<2.2$, and the effective
potential difference is 5%. Because $\Omega h^{2}$ are around $10^{-13}$, we
try with some value of $s$ and give the mass of new particles 5% difference,
the effective potential difference is also about 5% , it is found that $\Omega
h^{2}$ are almost unbiased.

We also make an extended comment, for looking at the strength of the EWPT when
gravitational waves are included. We have considered two parameters $S$ and
$\alpha$. A first-order EWPT is strong when $\alpha>1$ and $S>1$. So, in the
ZB model, the EWPT is the first-order one but not strong because $\alpha$ is
only about $10^{-2}$. In addition, we have further investigated the
gravitational wave spectrum over time in the phase transition by using the
formalism in Ref. grahamgw1 , but this requires knowing the exact decrease in
temperature over time or knowing the details of the expansion of the universe
and the growth of bubbles, which must not be approximate.

In this paper, we only calculate the contributions of components $m_{h^{\pm}}$
and $m_{k^{\pm\pm}}$. However the full masses of $h^{\pm}$ and $k^{\pm\pm}$
are Eq.(5), so that their masses still depend on $u_{1}$ and $u_{2}$.

To comply with some two-loop effects, for neutrino mass generation, the masses
of these new bosons are not allowed to exceed 2 TeV test . In addition to
$M_{k{\pm\pm}}>200$ GeV and $2M_{h^{\pm}}<M_{k^{\pm\pm}}$. With the results in
this paper, to have the first order phase transition and the sphaleron
energies are in the range of 5-10 TeV, their massed are about less than 300
GeV. However, this does not violate the above conditions because we have not
taken into account $u_{1}$ and $u_{2}$.

The recent new effects at the LHC have shown signals of di-bosons that have
masses around 2 TeV atlas . It has raised hopes for the existence of new
bosons in the ZB model. An interesting effect in this model is the small mass
generation of neutrino at the two-loop level. We think that this two-loop
effect will have a certain influence on the EWPT. Therefore, we will calculate
the effective potential at the two-loop to examine the baryogenesis problem
fully in the ZB model.

## Acknowledgement

In memory of the day VQP’s venerable teacher, Dr. Vo Thanh Van, passed away on
5 July 2016, this series of our papers highlights the continuation of what he
has done through many years in the Department of Theoretical Physics.

## References

  * (1) M.Fukugita and Yanagida, Phys. Lett. B 174, 45 (1986).
  * (2) A. D. Sakharov, JETP Lett.5, 24 (1967).
  * (3) M. I. Dine, R. G. Leigh, P. Huet, A. Linde, and D. Linde, Phys. Rev. D 46, 550 (1992).
  * (4) K. Kajantie, M. Laine, K. Rummukainen, and M. Shaposhnikov, Phys. Rev. Lett. 77, 2887 (1996); F. Csikor, Z. Fodor, and J. Heitger, Phys. Rev. Lett. 82, 21 (1999); J. Grant, M. Hindmarsh, Phys. Rev. D 64, 016002 (2001); M. D’Onofrio, K. Rummukainen, A. Tranberg, JHEP 08, 123 (2012).
  * (5) M. D’Onofrio, K. Rummukainen, A. Tranberg, Phys. Rev. Lett. 113, 141602 (2014).
  * (6) Vo Quoc Phong, Vo Thanh Van, and Hoang Ngoc Long, Phys. Rev. D 88, 096009 (2013), arXiv:1309.0355[hep-ph]; A. Menon, D. E. Morrissey, and C. E. M. Wagner, Phys. Rev. D 70, 035005 (2004); S. W. Ham, S. K. Oh, C. M. Kim, E. J. Yoo, and D. Son, Phys. Rev. D 70, 075001 (2004).
  * (7) M. Bastero-Gil, C. Hugonie, S. F. King, D. P. Roy, and S. Vempati, Phys. Lett. B 489, 359 (2000); A. Menon, D. E. Morrissey, and C. E. M. Wagner, Phys. Rev. D 70, 035005 (2004); S. W. Ham, S. K. Oh, C. M. Kim, E. J. Yoo, and D. Son, Phys. Rev. D 70, 075001 (2004).
  * (8) J. M. Cline, G. Laporte, H. Yamashita, S. Kraml, JHEP 0907, 040 (2009); Aleksandr Azatov, Miguel Vanvlasselaer, JHEP09, 085 (2020).
  * (9) S. Kanemura, Y. Okada, E. Senaha, Phys. Lett. B 606, 361-366 (2005); G. C. Dorsch, S. J. Huber, J. M. No, JHEP 10 (2013) 029.
  * (10) S. W. Ham, S-A Shim, and S. K. Oh, Phys. Rev. D 81, 055015 (2010).
  * (11) S. Das, P. J. Fox, A. Kumar, and N. Weiner, JHEP 1011, 108 (2010); D. Chung and A. J. Long, Phys. Rev. D 84, 103513 (2011); M. Carena, N. R. Shaha, and C. E. M. Wagner,Phys. Rev. D 85, 036003 (2012); A. Ahriche and S. Nasri, Phys. Rev. D 85, 093007 (2012); D. Borah and J. M. Cline, Phys. Rev. D 86, 055001 (2013).
  * (12) V. Q. Phong, H. N. Long, V. T. Van, N. C. Thanh, Phys. Rev. D 90, 085019 (2014), arXiv:1408.5657[hep-ph]; J. Sá Borges, R. O.Ramos, Eur. Phys. J. C 76: 344 (2016).
  * (13) V. Q. Phong, H. N. Long, V. T. Van, L. H. Minh, Eur. Phys. J. C 75,342 (2015), arXiv:1408.5657[hep-ph].
  * (14) J. R. Espinosa, T. Konstandin and F. Riva, Nucl. Phys. B 854, 592 (2012).
  * (15) S. Kanemura, E. Senaha, T. Shindou and T. Yamada, JHEP 1305, 066 (2013).
  * (16) D. J. H. Chung and A. J. Long, Phys. Rev. D 81, 123531 (2010).
  * (17) G. Barenboim and N. Rius, Phys. Rev. D 58, 065010 (1998).
  * (18) S. Profumo, M. J. Ramsey-Musolf, G. Shaughnessy, JHEP 0708 (2007) 010; S. Profumo, M. J. Ramsey-Musolf, C. L. Wainwright, P. Winslow, Phys. Rev. D 91, 035018 (2015); D. Curtin, P. Meade, C-T. Yu, JHEP 11, 127 (2014); M. Jiang, L. Bian, W. Huang, J. Shu, Phys. Rev. D 93, 065032 (2016).
  * (19) M. Carena, G. Nardini, M. Quiros, C. E.M. Wagner, Nucl. Phys. B 812:243-263 (2009); A. Katz, M. Perelstein, M. J. Ramsey-Musolf, P. Winslow, Phys. Rev. D 92, 095019 (2015); J. Kozaczuk, S. Profumo, L. S. Haskins, C. L. Wainwright, JHEP 1501, 144 (2015).
  * (20) H. H. Patel, M. J. Ramsey-Musolf, Phys.Rev. D. 88, 035013 (2012); N. Blinov, J. Kozaczuk, D. E. Morrissey, C. Tamarit, Phys. Rev. D 92, 035012 (2015); S. Inoue, G. Ovanesyan, M. J. Ramsey-Musolf, Phys. Rev. D 93, 015013 (2016).
  * (21) H. H. Patel, M. J. Ramsey-Musolf, JHEP 1107, 029 (2011); G. W. Anderson and L. J. Hall, Phys. Rev. D 45, 2685 (1992); H. H.Patel and M.J.Ramsey-Musolf, M. Garny and T.Konstandin, JHEP 1207, 189 (2012); C.W. Chiang, M.J. Ramsey-Musolf and E. Senaha, Phys. Rev. D 97, 015005 (2018).
  * (22) S. A.R. Ellis, S. Ipek, G. White, JHEP 1908, 002 (2019); J. De Vries, M. Postma, J. van de Vis, JHEP 1904, 024 (2019); J. de Vries, M. Postma, J. van de Vis , G. White, JHEP 1801, 089 (2018); C. Balazs, G. White, J. Yue, JHEP 1703, 030 (2017).
  * (23) A. Ahriche, Phys. Rev. D 75, 083522 (2007).
  * (24) A. Ahriche, Eur. Phys. J. C 66, 333 (2010); A. Ahriche, T. A. Chowdhury and S. Nasri, JHEP 1411, 096 (2014).
  * (25) A. Ahriche and S. Nasri, JCAP 1307, 035 (2013); A. Ahriche, G. Faisel, S. Y. Ho, S. Nasri and J. Tandean, Phys. Rev. D 92, 035020 (2015); A. Ahriche, K. L. McDonald and S. Nasri, Phys. Rev. D 92, 095020 (2015); A. Ahriche, S. M. Boucenna and S. Nasri, Phys. Rev. D 93, 075036 (2016); A. Ahriche, K. Hashino, S. Kanemura and S. Nasri, Phys. Lett. B 789, 119 (2019).
  * (26) K. Fuyuto and E. Senaha, Phys. Rev. D 90, 015015 (2014); K. Fuyuto and E. Senaha, Phys. Lett. B 747, 152 (2015); K. Funakubo and E. Senaha, Phys. Rev. D 79, 115024 (2009).
  * (27) M. Spannowsky and C. Tamarit, Phys. Rev. D 95, 015006 (2017).
  * (28) C. Grojean, G. Servant, J. D. Wells, Phys.Rev. D 71 036001 (2005).
  * (29) C. Delaunay, C. Grojean, J. D. Wells, JHEP 0804, 029 (2008).
  * (30) A. Kusenko, L. Pearce, and L. Yang, Phys. Rev. Lett. 114, 061302 (2015).
  * (31) V. Q. Phong, H. N. Long, V. T. Van, L. H. Minh, Eur. Phys. J. C 75, 342 (2015), arXiv:1409.0750[hep-ph].
  * (32) C. W. Chiang, T. Yamada, Phys. Let. B 735, 295 (2014).
  * (33) V. Q. Phong, N. C. Thao, H. N. Long, Phys. Rev. D 97, 115008 (2018), arXiv:1511.00579.
  * (34) J. R. Espinosa, T. Konstandin, J. M. No and M. Quiros, Phys. Rev. D 78, 123528 (2008).
  * (35) D. Comelli, J.R. Espinosa, Phys.Rev.D 55, 6253 (1997).
  * (36) M. Joyce, Phys. Rev. D 55, 1875 (1997).
  * (37) K. Funakubo and E. Senaha, Phys. Rev. D 79, 115024 (2009); M. Dvornikov and V. B. Semikoz, Phys. Rev. D 87, 025023 (2013); T. M. Gould and I. Z. Rothstein, Phys. Rev. D 48, 5917 (1993); K. Fuyuto and E. Senaha, Phys. Rev. D 90, 015015 (2014).
  * (38) S. Braibant, Y. Brihaye and J. Kunz, Int. J. Mod. Phys. A. 08, No. 31, 5563 (1993).
  * (39) K. Hagiwara, S. Ishihara, R. Szalapski and D. Zeppenfeld, Phys. Rev. D 48 2182 (1993).
  * (40) F. R. Klinkhamer and N. S. Manton, Phys. Rev. D 30, 2212 (1984).
  * (41) M. Spannowsky and C. Tamarit, Phys. Rev. D 95, 015006 (2017).
  * (42) S. Das, P. J. Fox, A. Kumar, and N. Weiner, JHEP 1011, 108 (2010).
  * (43) F. R. Klinkhamer and P. Nagel, Phys. Rev. D 96, 016006 (2017).
  * (44) A. Friedlander, I. Banta, J. M. Cline, and D. Tucker-Smith, Phys. Rev. D 103, 055020 (2021).
  * (45) S. W. Ham and S. K. Oh, Phys. Rev. D 70, 093007 (2004).
  * (46) X. Gan, A. J. Long, and L-T. Wang, Phys. Rev. D 96, 115018 (2017).
  * (47) Y. Brihaye and J. Kunz, Phys. Rev. D 48, 3884 (1993).
  * (48) C. L. Wainwright, Computer Physics Communications 183, 2006 (2012); C. L. Wainwright, Comput. Phys. Commun. 183, 2006 (2012).
  * (49) Vo Quoc Phong, Phan Hong Khiem, Ngo Phuc Duc Loc, Hoang Ngoc Long, Phys. Rev. D 101, 116010 (2020), arXiv:2003.09625[hep-ph].
  * (50) Kip S. Thorne, Gravitational Waves, arXiv:gr-qc/9506086v1 (1995).
  * (51) L. Leitao, A. Megevand, A. D. Sanchez, JCAP10, 024 (2012); M. S Tuner and S. Wilzek, Phy. Rev.65, 3080 (1990); C. Grojean, G. Servant, Phys.Rev.D75, 043507 (2007).
  * (52) S.V. Demidov, D. S. Gorbunov, D. V. Kirpichnikov, Physics Letters B 779, 191 (2018).
  * (53) S. J. Huber, T. Konstandin, JCAP 0809, 022 (2008); Ryusuke Jinno, Masahiro Takimoto, Phys. Rev. D 95, 024009 (2017).
  * (54) D. J. H. Chung, and P. Zhou, Phys. Rev. D 82, 024027 (2010).
  * (55) Zh. Kang, JHEP02, 115 (2018).
  * (56) R. Zhou, L. Bian and H-K Guo, Phys. Rev. D 101, 091903 (2020).
  * (57) B. Laurent and J. M. Cline, Phys. Rev. D 102, 063516 (2020).
  * (58) J. M. Cline, A. Friedlander, D-M He, K. Kainulainen, B. Laurent, D. Tucker-Smith, Phys. Rev. D 103, 123529 (2021).
  * (59) O. Gould, J. Kozaczuk, L. Niemi, M. J. Ramsey-Musolf, T. V. I. Tenkanen, D. J. Weir, Phys. Rev. D 100, 115024 (2019).
  * (60) M. F. Axen, S. Banagiri, A. Matas, C. Caprini, V. Mandic, Phys. Rev. D 98, 103508 (2018).
  * (61) M. Kamionkowski, A. KosowsKy and M. Turner,Phys. Rev. D 49, 2837 (1994).
  * (62) A. D. Dolgov, D. Grasso, A. Nicolis, Phys.Rev. D66, 103505 (2002).
  * (63) K. Schmitz, JHEP01, 097 (2021); R. Apredaa, M. Maggioreb, A. Nicolisc and A. Riotto, Nucl. Phys. B631, 342-368 (2002).
  * (64) C. Grojean and G. Servant, Phys. Rev. D75, 043507 (2007); M. Maggiore, Phys. Rept. 331, 283-367 (2000).
  * (65) A. Zee, Nucl. Phys. B 264, 99 (1986); K. S. Babu, Phys. Lett. B 203, 132 (1988).
  * (66) H. Okada, T. Toma, K. Yagyu, Phys. Rev. D 90, 095005 (2014); S. Baek, JHEP 08, 023 (2015); H. Hatanaka, K. Nishiwaki, H. Okada, Y. Orikasa, Nucl.Phys. B 894 (2015).
  * (67) G. W. Anderson and L. J. Hall, Phys. Rev. D 45, 2685-2698 (1992).
  * (68) Chiara Caprini, Mark Hindmarsh, Stephan Huber, Thomas Konstandin, Jonathan Kozaczuk, Germano Nardini, Jose Miguel No, Antoine Petiteau, Pedro Schwaller, Geraldine Servant, David J. Weir, JCAP04, 001 (2016).
  * (69) A. Kosowsky, M. S. Turner, and R. Watkins, Phys. Rev. D 45, 4514 (1992); A. Kosowsky, M. S. Turner, and R. Watkins, Phys. Rev. Lett. 69, 2026 (1992); A. Kosowsky and M. S. Turner, Phys. Rev. D 47, 4372 (1993); S. J. Huber and T. Konstandin, JCAP 08099, 022 (2008); R. Jinno and M. Takimoto, Phys. Rev. D 95, 024009 (2017); R. Jinno and M. Takimoto, JCAP 1901, 060 (2019).
  * (70) Huai-Ke Guo, K. Sinha, D. Vagie, G. White, JCAP 01, 001 (2021).
  * (71) T. Ohlsson, T. Schwetz, He Zhang, Phys. Lett. B 681, 269 (2009); Hoang Ngoc Long and Vo Van Vien, Int. J. Mod. Phys. A 29, 1450072 (2014).
  * (72) Huai-Ke Guo, K. Sinha, D. Vagie, G. White, JHEP06, 164 (2021); D. Croon, O. Gould, P. Schicho, T. V. I. Tenkanen, G. White, JHEP 04, 055 (2021).
  * (73) J. Herrero-Garcia, M. Nebot, N. Rius and A. Santamaria, Nucl. Phys. B 885, 542-570 (2014); J. Herrero-Garcia, M. Nebot, N. Rius, A. Santamaria, Proceedings of the 37th International Conference on High Energy Physics (ICHEP 2014).
  * (74) G. Aad et al, [ATLAS Collaboration], JHEP 12, 55 (2015); V. Khachatryan et al [CMS Collaboration], JHEP 1408, 173 (2014).


# Josephson effect in a fractal geometry

Morten Amundsen  Center for Quantum Spintronics, Department of Physics,
Norwegian University of Science and Technology, NO-7491 Trondheim, Norway
Vladimir Juričić  Departamento de Física, Universidad Técnica Federico Santa
María,
Casilla 110, Valparaíso, Chile  Nordita, KTH Royal Institute of Technology and
Stockholm University,
Hannes Alfvéns väg 12, SE-106 91 Stockholm, Sweden  Jabir Ali Ouassou
Department of Computer Science, Electrical Engineering and Mathematical
Sciences,
Western Norway University of Applied Sciences, NO-5528 Haugesund, Norway

###### Abstract

The Josephson effect is a hallmark signature of the superconducting state,
which, however, has been sparsely explored in non-crystalline superconducting
materials. Motivated by this, we consider a Josephson junction consisting of
two superconductors with a fractal metallic interlayer, which is patterned as
a _Sierpiński carpet_ by removing atomic sites in a self-similar and scale-
invariant manner. We here show that the fractal geometry has direct observable
consequences on the Josephson effect. In particular, we demonstrate that the
form of the supercurrent-magnetic field relation as the fractal generation
number increases can be directly related to the self-similar fractal geometry
of the normal metallic layer. Furthermore, the maxima of the corresponding
diffraction pattern directly encode the self-repeating fractal structure in
the course of fractal generation, implying that the corresponding magnetic
length directly probes the shortest length scale in the given fractal
generation. Our results should motivate future experimental efforts to verify
these predictions in designer quantum materials.

Fractals are paradigmatic non-crystalline structures featuring a non-integer
spatial dimension with a self-similarly repeating pattern at smaller and
smaller length scales, which are rather ubiquitous in Nature, from living
organisms to geological objects. Mandelbrot (1982) In the quantum realm, they
are realized most commonly in terms of wavefunctions exhibiting a fractal
behavior when electrons confined on a plane are subjected to a perpendicular
magnetic field, e.g. in the Hofstadter butterfly Hofstadter (1976), as well as
in disordered electronic systems. Soukoulis and Economou (1984); Schreiber and
Grussbach (1991) The interplay of fractal geometry in non-integer dimensions
and quantum-mechanical collective behavior of electrons living on fractals has
attracted renewed interest with advances in assembling of fractal structures
Shang  _et al._ (2015) and the experimental observation of the emerging
fractal features of the electronic wavefunctions. Kempkes  _et al._ (2019)
These observations spurred further interest in this problem, especially in
light of possible nontrivial topological properties of the electronic
wavefunctions on the fractals. Brzezińska, Cook, and Neupert (2018); Pai and
Prem (2019); Iliasov, Katsnelson, and Yuan (2020); Manna  _et al._ (2020);
Fremling  _et al._ (2020); Yang  _et al._ (2020); Manna  _et al._ (2022);
Manna, Nandy, and Roy (2022); Ivaki  _et al._ (2022); Manna and Roy (2023)
Quite surprisingly, however, superconductivity in this context has been rather
sparsely explored. Hori  _et al._ (2024); Iliasov, Katsnelson, and Bagrov
(2024); Meyer  _et al._ (2002); Kingni  _et al._ (2019); Manna, Das, and Roy
(2022) Intriguingly, a very recent study has however suggested that artificial
fractal geometries can enhance the critical temperature $T_{c}$ of
superconductors.Iliasov, Katsnelson, and Bagrov (2024)

Superconductivity is a macroscopic quantum phenomenon, where a large number of
electron pairs all condense into the same macroscopic quantum state. This
quantum state can be described using a macroscopic wave function
$\Psi(\bm{r},t)=|\Psi(\bm{r},t)|e^{i\varphi(r,t)}$. When two different
superconductors are connected to the same metal, their wave functions can leak
into the metal and interfere with each other. This can result in the
_Josephson effect_ ,Josephson (1962); Golubov, Kupriyanov, and Il’ichev (2004)
whereby a current that depends on the phase difference $\delta\varphi$ between
the wavefunctions in the two superconductors flows through the normal metal.
Josephson junctions have found numerous applications, ranging from precise
metrology (e.g. the SQUID) to novel computing paradigms (e.g. RSFQ logic or
phase qubits).

Figure 1:  Geometry of the S/N/S Josephson junction with the N layer
exhibiting a fractal (self-similar) structure. The setup consists of two
superconductors (yellow) connected by a normal metal (black) that has been
patterned as a Sierpiński fractal by recursively creating insulating holes
according to eq. 3. A homogeneous in-plane magnetic field
$\bm{B}\sim\bm{e}_{x}$ is then applied to the Josephson junction, causing
spatially varying screening currents $\bm{J}\sim\bm{e}_{z}$ to flow between
the superconductors. Depending on the magnitude $B$ of the applied magnetic
field, the screening currents in different parts of the interlayer may either
interfere constructively or destructively, causing a kind of interference
pattern to appear in the current–field relation $I(B)$.

The Josephson effect in fractal geometries has received only very limited
attention in the literature so far. Meyer  _et al._ (2002); Kingni  _et al._
(2019); Kruchinin  _et al._ (2005) In this paper, we fill this gap by
considering a heterostructure made of a sandwich of conventional $s$-wave
superconductors (S) and a normal material (N) with the Sierpiński fractal
geometry. Note that in contrast to e.g. Ref. Meyer  _et al._ , 2002, we do not
here consider fractal properties at the level of Josephson junction _arrays_ ,
but rather a fractal structure within _one_ junction. Moreover, while e.g.
Refs. Kingni  _et al._ , 2019; Kruchinin  _et al._ , 2005 consider interlayers
with fractal characteristics that are naturally acquired during lithography,
we here consider artificial fractal geometries that can exhibit much more
dramatic phenomenology. In particular, we show that fractal geometry, here
taken to be a Sierpiński carpet (fig. 1), has direct observable consequences
on the Josephson effect. The form of the supercurrent-magnetic field relation
in the process of the fractal generation, as we show, can be directly related
to the self-similar fractal pattern of the normal metallic layer, as shown in
fig. 2. Furthermore, the corresponding diffraction pattern, through its
maxima, encodes the self-repeating fractal structure in the course of fractal
generation (fig. 3), which implies that the corresponding magnetic length
directly probes the shortest length scale in the given fractal generation. Our
findings should motivate experimental pursuits to probe our theoretical
predictions in the designer quantum materials.

We consider an S/N/S Josephson junction as illustrated in fig. 1.Josephson
(1962) We take the stack to be oriented along the $z$ axis, and to have
dimensions $W\times W$ in the $xy$ plane. Each S layer is a conventional BCS
superconductor,Bardeen, Cooper, and Schrieffer (1957) and is assumed to be
much longer than the magnetic penetration depth $\lambda$ along the $z$
direction. The N layer is a thin layer of a nominally normal metal, but has in
our case been patterned as a Sierpiński carpet: It contains a fractal pattern
of holes, where atoms have been removed from the otherwise normal metal in a
self-similar and scale-invariant manner (see fig. 1). In an experiment, it is
likely easier to deposit a thin layer of an electric insulator at the
locations of the holes in our geometry, which would have the same effect as
our missing lattice sites. Artificial fractal geometries in condensed matter
systems—including Sierpiński fractals—have in recent years been experimentally
realized.Kempkes  _et al._ (2019)

The Josephson effect between the two superconductors produces dissipationless
currents that flow through the fractal layer whenever a phase difference is
enforced. We assume that the Josephson penetration depth is much larger than
the in-plane dimensions of the system, and hence any self-screening effects
may be neglected, meaning that the currents flow solely in the $z$ direction.
The resulting current density therefore takes the form

$\bm{J}(x,y)=J(x,y)\,\bm{e}_{z}=J_{0}g(x,y)\sin[\varphi(x,y)]\,\bm{e}_{z},$
(1)

where $\delta\varphi$ is the gauge-invariant phase difference between two
vertically separated points deep inside each superconductor, and $g(x,y)$ is a
form factor, indicating the shape of the fractal. To leading order, we have
here assumed that the critical current density $J(x,y)$ is uniform and equal
to $J_{0}$ at coordinates without a hole, whereas it drops to zero in the
insulating holes. The critical current density is thus given as
$J_{c}(x,y)=J_{0}g(x,y)$. If we express the coordinates $x=ia$ and $y=ja$ in
terms of a lattice constant $a$, we can then write

$g(x,y)=S^{N}_{ij},$ (2)

where $S^{N}_{ij}$ are matrix elements of an $N$’th-order Sierpiński carpet.
The Sierpiński carpet is a fractal structure that has the properties of self-
similarity and scale invariance—at least down to the scale of individual unit
cells—and is illustrated in fig. 1. Such a lattice may be constructed using
the recursive algorithm

$\displaystyle\bm{S}^{N}$ $\displaystyle=\bm{S}^{N-1}\otimes\bm{S}^{1},$
$\displaystyle\bm{S}^{1}$ $\displaystyle=\begin{pmatrix}\,1&1&1\;\\\
\,1&0&1\;\\\ \,1&1&1\;\end{pmatrix},$ $\displaystyle\bm{S}^{0}$
$\displaystyle=1,$ (3)

where $\otimes$ is the Kronecker matrix product. This implies that we can also
write the $N^{\rm th}$ order Sierpiński carpet
$\bm{S}^{N}=(\bm{S}^{1})^{\otimes N}$ compactly as a Kronecker matrix power.

Figure 2: Fourier transform $J_{c}(k_{x},k_{y})$ of the current distribution
in a Sierpiński Josephson junction (upper panel). We here considered a 6th-
order Sierpiński lattice as the interlayer. According to the Dynes–Fulton
method, the net current $I(\beta)$ that arises for an applied magnetic field
$B\sim\beta$ along the $x$ axis is found from the Fourier-transformed current
density along the $y$ axis, as indicated by the blue line in the upper panel,
with the result shown in the lower panel.

We now consider a magnetic field $\bm{B}=B\bm{e}_{x}$ applied in the thin-film
plane of the fractal interlayer. This produces screening currents in the $yz$
plane of the superconductors, within a distance on the order of their
penetration depth $\lambda$ away from the interface to the normal metal. These
currents take the form of elongated loops that flow purely vertically through
the thin normal metal, and give rise to a position-dependence in the gauge-
invariant phase difference, $\phi=\phi_{0}-\frac{2e}{\hbar}Bdy$, where
$\phi_{0}$ is an externally enforced phase difference, and $d=l+2\lambda$ is
the effective thickness of the junction in the $z$ direction, with $l$ the
thickness of the normal metal. The result are spatially varying current
contributions in the normal metal that can interfere constructively or
destructively depending on the precise value of the magnetic field $B$. The
total current flowing through the system at a given $B$ is given as the
surface integral over the cross section, which is quadratic with widths $W$,

$I(B,\phi_{0})=J_{0}\int_{-W/2}^{W/2}\int_{-W/2}^{W/2}dxdy\;g(x,y)\sin\left[\phi_{0}-\frac{2e}{\hbar}Bdy\right],$
(4)

As was first observed by Dynes and Fulton,Dynes and Fulton (1971) this
expression is simply the imaginary part of the Fourier transform of the
current density, which in our case becomes the Fourier transform of the
fractal form factor, $g(k_{x},k_{y})$. Furthermore, since $g(x,y)$ has
inversion symmetry, $g(k_{x},k_{y})$ is purely real. The critical current is
thus found at $\phi=\pi/2$, and is given as

$I_{c}(B)=|J_{c}(k_{x}=0,k_{y}=\beta)|=J_{0}\int d\bm{r}\;g(\bm{r})e^{-i\beta
y},$ (5)

with $\bm{r}=x\bm{e}_{x}+y\bm{e}_{y}$. Here, the magnetic field $B$ is
parametrized via $\beta=(2\pi/W)(\Phi/\Phi_{0})$, where $\Phi=BWd$ is the net
flux passing through the central parts of the Josephson junction, and
$\Phi_{0}=h/2e$ is the flux quantum.

Using the approach outlined above, the current–field relation $I_{c}(B)$ can
be extracted as follows. First, an $N^{\rm th}$ order Sierpiński carpet is
constructed recursively using eq. 3. This is subsequently employed to find
$J_{c}(x,y)$ according to eqs. 1 and 2, which can then be run through a 2D
Fast Fourier Transform (FFT) algorithm to obtain $J_{c}(k_{x},k_{y})$.
Finally, we can then simply extract the critical current along a line in the
2D plane, as determined by the direction of the applied field and the chosen
gauge. This procedure, and the resulting $I_{c}(B)$, are shown in fig. 2 for a
junction with a 6th-order Sierpiński interlayer. Unsurprisingly, the Fourier
transform of a fractal contains self-similar features also, and ends up having
structures on all length scales from $\pi/W$ to $\pi/a$. The resulting field
dependence has a very different structure from the well-known Fraunhofer and
SQUID patterns that arise in comparable non-fractal junctions.

Figure 3: The critical current $I_{c}$ as a function of magnetic field,
$N_{\Phi}=BdW/\Phi_{0}$, for various generations $N$ of the Sierpiński carpet.
(a) Line plots for $N\in\left\\{3,4,5\right\\}$, shifted vertically for
clarity. The vertical lines indicate the location of new peaks introduced by
an increase in $N$. (b) A surface plot of $I_{c}(B)$ for a wider range of
values of $N$ and $N_{\Phi}$. The dashed line is a guide for the eye, showing
approximately where $I_{c}(B)$ becomes independent of $N$.

We can arrive at this result by purely analytical means. This can be done by
treating current as a parallel coupling of the current passing through each of
the constituent elements of a generation of the Sierpiński carpet. To
illustrate this we consider, for instance, the first generation, as is
obtained from fig. 1 by retaining only the largest hole in the center. In that
case, the total current may be split into three parts, two of which form
rectangular junctions with dimensions $W\times W/3$, and one part which has
the form of a SQUID formed by two junctions with dimensions $W/3\times W/3$
and an equally sized hole in between. The rectangular contributions provide a
current in a Fraunhofer pattern, $J_{0}\frac{W^{2}}{3}\mathrm{sinc}\left(\pi
N_{\Phi}\right)$, whereas the SQUID contribution gives a current
$2J_{0}\left(\frac{W}{3}\right)^{2}\mathrm{sinc}\left(\frac{\pi
N_{\Phi}}{3}\right)\cos\left(\frac{2\pi N_{\Phi}}{3}\right)$. Here,
$N_{\Phi}=BdW/\Phi_{0}$ is the number of flux quanta passing through the
system, and $\mathrm{sinc}(x)=\sin(x)/x$. Adding up each contribution gives
the critical current for the first-generation Sierpiński carpet

$I^{(1)}_{c}(B)=\frac{J_{0}A_{F}^{(1)}}{4}\left|\mathrm{sinc}\left(\frac{\pi
N_{\Phi}}{3}\right)\left[1+3\cos\frac{2\pi N_{\Phi}}{3}\right]\right|,$

where $A_{F}^{(N)}=\left(\frac{8}{9}\right)^{n}W^{2}$ is the surface area of
the fractal at generation $n$. An expression for arbitrary generation $N$ can
be obtained by repeating this process recursively, replacing each uniform sub-
square with the shape of the previous generation. This leads to

$I^{(N)}_{c}(B)=\frac{J_{0}A^{(N)}_{F}}{4^{N}}\left|\,\mathrm{sinc}\left(\frac{\pi
N_{\Phi}}{3^{N}}\right)\prod_{n=1}^{N}\left[1+3\cos\left(\frac{2\pi
N_{\Phi}}{3^{n}}\right)\right]\,\right|,$ (6)

The 0th-order Sierpiński carpet is equivalent to a uniform square lattice, and
indeed we regain the conventional Fraunhofer pattern for $N=0$ when no terms
appear in the right-hand product, $I^{(0)}_{c}(B)=J_{0}W^{2}\mathrm{sinc}(\pi
N_{\Phi})$. In the opposite limit $N\to\infty$ of a high-order Sierpiński
carpet, the prefactor approaches $\mathrm{sinc}(0)=1$. The current–field
relation then reduces to a product of similar factors $1+3\cos(\beta W_{n})$
arising from the Sierpiński pattern on different scales $W_{n}=W/3^{n}$. This
result shows clearly that $I_{c}(B)$ exhibits self-similarity across different
length scales, which arises from the fractal geometry of the Josephson
junction interlayer. This self-similarity can be seen in the 2D Fourier
transform in fig. 2, where it manifests as superimposed lattices of
$3^{n}\times 3^{n}$ squares.

The presence of multiple length scales has additional nontrivial observable
consequences as, for instance, the sensitivity of the diffraction pattern on
the order of fractal generation. In fig. 3(a) the diffraction pattern is shown
for a selection of fractal generations, $N\in\left\\{3,4,5\right\\}$. It can
be seen that each new generation introduces new features in the diffraction
pattern which were not present in the previous one. These features appear in
limited ranges of $N_{\Phi}$, and are characterized by a single peak of
constant height, the location of which is indicated by the vertical lines in
fig. 3(a). Indeed, for a fixed $N$, we observe that whenever $N_{\Phi}=3^{q}$
for an integer $q$, the square brackets in eq. 6 become constant and maximal
for all $n\leq q$. For $q<N$, $I_{c}^{(N)}$ increases with $q$, whereas for
$q\geq N$ it vanishes. The result is a local maximum at $q=N-1$, at which
point the critical current becomes
$I_{c}^{(N)}=3\sqrt{3}J_{0}A_{F}^{(N)}/16\pi$. The peak occurring at
$N_{\Phi}=3^{N-1}$ is a manifestation of the fact that the magnetic length
scale becomes small enough to resolve the lowest length scales of the fractal.
Alternatively, one may say that when fractal generation is increased from $N$
to $N+1$, the smallest length scale is reduced by a factor of 3. Hence, a flux
density corresponding to the rescaled magnetic field, $3B$, in the latter
(generation $N+1$), produces a similar response as $B$ in the former case
(generation $N$). Another interesting and related feature is that the critical
current for any $N_{\Phi}$ eventually becomes independent of $N$. This is
demonstrated in fig. 3(b), which shows a surface plot of the critical current
as a function of $N$ and $N_{\Phi}$. The dashed line indicates approximately
where this saturation occurs. For $N=7$, the diffraction pattern remains
constant at least for $N_{\Phi}\leq 1000$. This behavior is once again a
result of the phenomenon that as $N$ increases, smaller and smaller length
scales are introduced, which eventually become smaller than the magnetic
length scale. Their features cannot be resolved by the applied flux, and so
the system remains uninfluenced.

Here, we have focused on the specific case where the magnetic field
$\bm{B}\sim B\bm{e}_{x}$ is applied along the $x$ axis. Note however that a
similar Dynes–Fulton analysis is applicable for a magnetic field
$\bm{B}(\theta)=B(\cos\theta\,\bm{e}_{x}+\sin\theta\,\bm{e}_{y})$ oriented
along any direction in the $x-y$ plane, where each such direction would
provide a slightly different current–field relation $I_{c}(B\,|\,\theta)$. In
principle, one could in an experiment use the $I_{c}(B\,|\,\theta)$ curves
obtained for different field directions $\theta$ to reconstruct the 2D Fourier
transform $J_{c}(k_{x},k_{y})$ of the current density, and thus to obtain the
exact current distribution $J_{c}(x,y)$ in the interlayer. Such an approach
might be particularly interesting in Josephson junctions where the interlayer
might not be an artificially constructed fractal as considered here, but
rather a naturally occurring fractal or quasicrystal. In that case, this kind
of Dynes–Fulton analysis might provide a way to validate and elucidate the
non-crystalline, fractal or quasicrystalline, structure of the interlayer.

In this paper, we have considered a fractal Josephson junction composed of two
superconductors separated by a Sierpiński-patterned normal metal. We have
analytically calculated the current–field relation $I_{c}(B)$ that arises in
such junctions, and shown that this experimental signature retains a crucial
property of fractals: self-similarity across length scales. In particular, the
obtained diffraction pattern is directly related to the self-similar fractal
geometry probed through the magnetic length. Our results should be
consequential for experiments, particularly in the designer quantum materials,
where the realization of the normal fractal metallic layer should be feasible.

M.A. was supported by the Research Council of Norway through its Centres of
Excellence funding scheme Grant No. 262633 “QuSpin”. V.J. acknowledges the
support of the Swedish Research Council Grant No. VR 2019-04735 and Fondecyt
(Chile) Grant No. 1230933.
Data Availability
Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

## References

  * Mandelbrot (1982) B. B. Mandelbrot,  _The fractal geometry of Nature_ (W.H. Freeman and Co., 1982).
  * Hofstadter (1976) D. R. Hofstadter, “Energy levels and wave functions of bloch electrons in rational and irrational magnetic fields,” Phys. Rev. B 14, 2239–2249 (1976).
  * Soukoulis and Economou (1984) C. M. Soukoulis and E. N. Economou, “Fractal character of eigenstates in disordered systems,” Physical Review Letters 52, 565–568 (1984).
  * Schreiber and Grussbach (1991) M. Schreiber and H. Grussbach, “Multifractal wave functions at the Anderson transition,” Physical Review Letters 67, 607–610 (1991).
  * Shang  _et al._ (2015) J. Shang, Y. Wang, M. Chen, J. Dai, X. Zhou, J. Kuttner, G. Hilt, X. Shao, J. M. Gottfried, and K. Wu, “Assembling molecular Sierpiński triangle fractals,” Nature Chemistry 7, 389–393 (2015).
  * Kempkes  _et al._ (2019) S. N. Kempkes, M. R. Slot, S. E. Freeney, S. J. M. Zevenhuizen, D. Vanmaekelbergh, I. Swart, and C. M. Smith, “Design and characterization of electrons in a fractal geometry,” Nature Physics 15, 127–131 (2019).
  * Brzezińska, Cook, and Neupert (2018) M. Brzezińska, A. M. Cook, and T. Neupert, “Topology in the Sierpiński–Hofstadter problem,” Physical Review B 98, 205116 (2018).
  * Pai and Prem (2019) S. Pai and A. Prem, “Topological states on fractal lattices,” Physical Review B 100, 155135 (2019).
  * Iliasov, Katsnelson, and Yuan (2020) A. A. Iliasov, M. I. Katsnelson, and S. Yuan, “Hall conductivity of a Sierpiński carpet,” Physical Review B 101, 045413 (2020).
  * Manna  _et al._ (2020) S. Manna, B. Pal, W. Wang, and A. E. B. Nielsen, “Anyons and fractional quantum Hall effect in fractal dimensions,” Physical Review Research 2, 023401 (2020).
  * Fremling  _et al._ (2020) M. Fremling, M. van Hooft, C. M. Smith, and L. Fritz, “Existence of robust edge currents in sierpiński fractals,” Physical Review Research 2, 013044 (2020).
  * Yang  _et al._ (2020) Z. Yang, E. Lustig, Y. Lumer, and M. Segev, “Photonic floquet topological insulators in a fractal lattice,” Light: Science & Applications 9, 128 (2020).
  * Manna  _et al._ (2022) S. Manna, C. W. Duncan, C. A. Weidner, J. F. Sherson, and A. E. B. Nielsen, “Anyon braiding on a fractal lattice with a local Hamiltonian,” Physical Review A 105, L021302 (2022).
  * Manna, Nandy, and Roy (2022) S. Manna, S. Nandy, and B. Roy, “Higher-order topological phases on fractal lattices,” Physical Review B 105, L201301 (2022).
  * Ivaki  _et al._ (2022) M. N. Ivaki, I. Sahlberg, K. Pöyhönen, and T. Ojanen, “Topological random fractals,” Communications Physics 5, 327 (2022).
  * Manna and Roy (2023) S. Manna and B. Roy, “Inner skin effects on non-Hermitian topological fractals,” Communications Physics 6, 10 (2023).
  * Hori  _et al._ (2024) M. Hori, T. Sugimoto, T. Tohyama, and K. Tanaka, “Self-consistent study of topological superconductivity in two-dimensional quasicrystals,”  (2024), arXiv:2401.06355 .
  * Iliasov, Katsnelson, and Bagrov (2024) A. A. Iliasov, M. I. Katsnelson, and A. A. Bagrov, “Strong enhancement of superconductivity on finitely ramified fractal lattices,”  (2024), arXiv:2310.11497 .
  * Meyer  _et al._ (2002) R. Meyer, S. E. Korshunov, C. Leemann, and P. Martinoli, “Dimensional crossover and hidden incommensurability in Josephson junction arrays of periodically repeated Sierpiński gaskets,” Physical Review B 66, 104503 (2002).
  * Kingni  _et al._ (2019) S. T. Kingni, G. F. Kuiate, V. K. Tamba, A. V. Monwanou, and J. B. C. Orou, “Analysis of a fractal Josephson junction with unharmonic current–phase relation,” Journal of Superconductivity and Novel Magnetism 32, 2295–2301 (2019).
  * Manna, Das, and Roy (2022) S. Manna, S. K. Das, and B. Roy, “Noncrystalline topological superconductors,” (2022), arXiv:2207.02203 .
  * Josephson (1962) B. Josephson, “Possible new effects in superconductive tunnelling,” Physics Letters 1, 251–253 (1962).
  * Golubov, Kupriyanov, and Il’ichev (2004) A. A. Golubov, M. Y. Kupriyanov, and E. Il’ichev, “The current–phase relation in Josephson junctions,” Reviews of Modern Physics 76, 411–469 (2004).
  * Kruchinin  _et al._ (2005) S. P. Kruchinin, V. F. Klepikov, V. E. Novikov, and D. S. Kruchinin, “Nonlinear current oscillations in the fractal Josephson junction,” Materials Science Poland 23, 1009 (2005).
  * Bardeen, Cooper, and Schrieffer (1957) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, “Theory of superconductivity,” Physical Review 108, 1175–1204 (1957).
  * Dynes and Fulton (1971) R. C. Dynes and T. A. Fulton, “Supercurrent density distribution in Josephson junctions,” Physical Review B 3, 3015–3023 (1971).