# The N2HDM, Entropy Production and Stochastic Gravitational Waves

Arnab Chaudhuria***e-mail<EMAIL_ADDRESS>and Kazunori
Kohria,b,c,d†††e-mail<EMAIL_ADDRESS>
a Division of Science, National Astronomical Observatory of Japan,
Mitaka, Tokyo 181-8588, Japan.
b The Graduate University for Advanced Studies (SOKENDAI), Mitaka, Tokyo
181-8588, Japan.
c Theory Center, IPNS, and QUP (WPI), KEK, 1-1 Oho, Tsukuba, Ibaraki 305-0801,
Japan
d Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583,
Japan

###### Abstract

This study undertakes a reconsideration of the potential for a first-order
electroweak phase transition, focusing on the next-to-minimal two Higgs
doublet model (N2HDM). Our exploration spans diverse parameter spaces
associated with the phase transition, with a particular emphasis on examining
the generation of stochastic Gravitational Waves (GW) resulting from this
transition. The obtained results are meticulously compared against data from
prominent gravitational wave observatories, and the possibility of their
detection in the future GW observations have been established. In passing by
we analyse the strength of the phase transition through the production of
entropy during the electroweak phase transition.

Keywords: N2HDM, Entropy production, Gravitational Waves, Electroweak Phase
transition.

##  1 Introduction

While the Standard Model (SM) of particle physics has achieved remarkable
success in elucidating three of the four fundamental forces—namely,
electromagnetic, weak, and strong forces—it remains incomplete in its quest to
comprehensively describe all known elementary particles and the entirety of
fundamental interactions. Despite its popularity and widespread acceptance,
the SM grapples with certain unresolved phenomena. Notably, it fails to
provide a complete explanation for baryon asymmetry, lacks a cohesive
incorporation of the gravitational theory outlined in general relativity, and
does not address the universe’s accelerating expansion, a phenomenon possibly
governed by dark energy. Furthermore, the SM lacks a suitable candidate for
dark matter, leaving crucial aspects of our understanding of the cosmos
unaccounted for.

The origin of the observed asymmetry between matter and antimatter, known as
baryogenesis, adheres to Sakharov’s principle [1], encompassing three
essential conditions: i) Non-conservation of Baryon numbers, ii) Breaking of C
and CP invariance, and iii) Deviation from thermal equilibrium. To elucidate
the baryon asymmetry in the universe through baryogenesis, a pivotal
prerequisite is a strong first-order electroweak phase transition (EWPT)
during the early universe. The cosmic EWPT transpired as the primeval universe
cooled sufficiently, causing the Higgs field’s potential to settle at a non-
zero minimum. Consequently, the symmetry of the theory $SU(3)_{C}\times$
$SU(2)_{L}\times U(1)_{Y}$ broke into $U(1)_{\rm em}$. At the onset of this
first-order EWPT, bubbles of the broken phase emerged, generating baryon-
antibaryon asymmetry outside the bubble walls.

However, with the subsequent discovery of the SM Higgs boson [2], it became
apparent that the EWPT within the SM framework is a smooth crossover phase
transition (see [3, 4]). Consequently, for a successful electroweak
baryogenesis (EWBG), a first-order phase transition theory beyond the SM is
imperative. Studies have delved into strong first-order phase transitions in
$Z_{1}$, $Z_{2}$ and $Z_{3}$ extensions of the SM [5, 6, 7, 8, 9] to explore
viable alternatives.

In essence, the strength of the EWPT hinges on the intricate interplay between
the scalar potential’s behavior at high and low temperatures. The critical
temperature, determined through computation, serves as a pivotal parameter,
revealing shifts in the global minimum within the scalar potential concerning
the temperature (T) of the Universe. To gain a precise understanding of the
electroweak phase transition (EWPT), it becomes imperative to delve into the
dynamics of bubble nucleation. This process, essential for the progression of
the phase transition, involves the formation of bubbles [10]. Particularly
during the first-order EWPT, the dynamics of bubble nucleation can give rise
to stochastic gravitational waves (GWs) in the early Universe [11, 12, 13, 14,
15, 16]. The detectability of these GWs becomes a focal point, as they hold
the potential to be observed in various gravitational wave experiments. Not
only the production of stochastic GW is a pivotal point for a first order EWPT
but also the high production/release of entropy into the plasma which can
dilute any preexisting frozen out species like dark matter, baryon asymmetry,
etc. is a matter of severe importance. This release of entropy is almost
negligible within the framework of the SM because EWPT is a of crossover in
nature, [3]. But once the SM is extended, even by a singlet, there is a
significant production of entropy during EWPT, [5, 17].

This pursuit of gravitational wave detection has been a longstanding practice,
especially in the context of exploring diverse beyond the Standard Model (BSM)
frameworks [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. Notably, it serves as
a complementary avenue alongside traditional collider searches [29, 30, 31],
offering a multifaceted approach to unravel the mysteries of theoretical
frameworks beyond the Standard Model.

In this work we consider the next-to-minimal two Higgs doublet model (N2HDM),
which is a singlet extension to the two Higgs doublet model (2HDM), as the
extension to the SM. By scanning the parameter space for the first order phase
transition (FOPT) at the electroweak scale we study the production of the
stochastic GW and predict their possible detection in the future GW
observations.

This paper is organised as follows. In Section 2, we explain the scalar sector
of the N2HDM potential. Then, in Section 3, we discuss the thermal and one-
loop corrections to the scalar potential. Section 4 is divided into two
subsections. Subsection 4.1 deals with the entropy production scenarios due to
EWPT within the framework of the N2HDM and subsection 4.2 focuses on the study
of the stochastic gravitational waves production due to EWPT. This section is
followed by the conclusion 5, we summarize our analysis.

##  2 Model

Extensions to the Higgs sector, such as the two-Higgs-doublet model (2HDM) or
the Next-to-Minimal-2HDM (N2HDM), offer a diverse collider phenomenology,
deeply intertwined with early Universe physics [32, 33, 34, 35, 36, 37, 38,
39, 40, 41, 42, 43]. In this study, we specifically investigate this
connection within the framework of the N2HDM, with a focus on type II. Our
analysis delves into the thermal history, elucidating the evolution of Higgs
fields in the early Universe. We reveal that across significant portions of
the N2HDM parameter space, the thermal history diverges notably from the
conventional scenario of electroweak symmetry breaking occurring at an early
Universe temperature $T$ of around $\mathcal{O}(100)~{}\rm{GeV}$. Two critical
phenomena, namely electroweak symmetry non-restoration and vacuum trapping,
emerge as pivotal factors shaping this distinct thermal evolution.

The N2HDM tree-level scalar potential consists of two $\text{SU(2)}_{L}$ Higgs
doublets $\Phi_{1}$ and $\Phi_{2}$ and the real singlet field $\Phi_{S}$. A
$\mathbb{Z}_{2}$ symmetry is imposed, which is explicitly broken and
transforms the fields as
$\Phi_{1}\rightarrow\Phi_{1},\,\Phi_{2}\rightarrow-\Phi_{2},\,\Phi_{S}\rightarrow\Phi_{S}$.
At this moment we do not impose the extra $\mathbb{Z}^{\prime}_{2}$ symmetry
as mentioned in [36, 44]. This allows us to retain the cubic terms involving
$\Phi_{S}$ in the scalar potential, which takes the form:

$\displaystyle V_{\text{tree}}$
$\displaystyle=m_{11}^{2}\left|\Phi_{1}\right|^{2}+m_{22}^{2}\left|\Phi_{2}\right|^{2}-m_{12}^{2}\left(\Phi_{1}^{\dagger}\Phi_{2}+\text{h.c.}\right)+\frac{\lambda_{1}}{2}\left(\Phi_{1}^{\dagger}\Phi_{1}\right)^{2}+\frac{\lambda_{2}}{2}\left(\Phi_{2}^{\dagger}\Phi_{2}\right)^{2}$
$\displaystyle+\lambda_{3}\left(\Phi_{1}^{\dagger}\Phi_{1}\right)\left(\Phi_{2}^{\dagger}\Phi_{2}\right)+\lambda_{4}\left(\Phi_{1}^{\dagger}\Phi_{2}\right)\left(\Phi_{2}^{\dagger}\Phi_{1}\right)+\frac{\lambda_{5}}{2}\left[\left(\Phi_{1}^{\dagger}\Phi_{2}\right)^{2}+\mathrm{h.c.}\right]$
$\displaystyle+\frac{1}{2}m_{S}^{2}\Phi_{S}^{2}+a_{1}^{3}\Phi_{S}-a_{2}\Phi_{S}^{3}+a_{3}\left(\Phi_{1}^{\dagger}\Phi_{1}\right)\Phi_{S}$
$\displaystyle+\frac{\lambda_{6}}{8}\Phi_{S}^{4}+\frac{\lambda_{7}}{2}\left(\Phi_{1}^{\dagger}\Phi_{1}\right)\Phi_{S}^{2}+a_{4}\left(\Phi_{2}^{\dagger}\Phi_{2}\right)\Phi_{S}+\frac{\lambda_{8}}{2}\left(\Phi_{2}^{\dagger}\Phi_{2}\right)\Phi_{S}^{2}.$
(1)

It is clear from (2) that $a_{i}~{}(i=1,2,3,4)$ has the dimension of
mass/temperature. The $\mathbb{Z}_{2}$ symmetry in (2), serves to safeguard
against flavor-changing neutral currents (FCNCs) at the lowest order when
extended to the Yukawa sector. However, this symmetry is softly broken by the
presence of the $m_{12}^{2}$ term. In the third and fourth line of the tree-
level potential (2), we incorporate the impact of the singlet field. We
explore the scenario where $\Phi_{S}$ obtains a vacuum expectation value
(vev). To analyze this further, we expand the fields around the EW minimum as
follows:

$\Phi_{1}=\left(\begin{array}[]{c}\phi_{1}^{+}\\\
\frac{1}{\sqrt{2}}\left(v_{1}+\rho_{1}+i\eta_{1}\right)\end{array}\right),\quad\Phi_{2}=\left(\begin{array}[]{c}\phi_{2}^{+}\\\
\frac{1}{\sqrt{2}}\left(v_{2}+\rho_{2}+i\eta_{2}\right)\end{array}\right),\quad\Phi_{S}=v_{S}+\rho_{3},$
(2)

where $v_{1}$, $v_{2}$ and $v_{S}$ are the field vevs for the Higgs doublets
and the singlet field, respectively, at zero temperature. The fields
$\phi_{i}^{+}$ (where $i=1,2$) represent the complex scalar components linked
to the two Higgs doublets. The CP-even and odd fields are denoted by
$\rho_{I}$ (with $I=1,2,S$) and $\eta_{i}$, respectively. The vevs of the
doublets are represented by $v_{1}$ and $v_{2}$, defining the electroweak (EW)
scale $v=\sqrt{v_{1}^{2}+v_{2}^{2}}$, which is approximately $246$ GeV.

With the introduction of another discrete $\mathbb{Z}^{\prime}_{2}$ symmetry,
the $a_{i}$ vanishes and is the best “coordinate frame” for the singlet as it
makes explicitly apparent such symmetry. In case the singlet does not acquire
a vev, then it can be a suitable candidate for a dark matter particle but we
do not divulge into the study of dark matter.

At this point, it is important to lay out some differences between the 2HDM
and the N2HDM. Not only does the N2HDM has more degrees of freedom compared to
the 2HDM because of the extra singlet scalar but it exhibits a different
vacuum structure as well. At tree level, the general vacuum structure of the
2HDM allows for three different possible vacua that are given by the normal
EW-breaking vacuum, a CP-breaking and a charge-breaking (CB) vacuum [45, 46].
Given that this statement might not remain valid at higher orders or could be
disrupted by finite temperature effects, we entertain the possibility of a
more versatile vacuum structure, a characteristic inherent to the N2HDM.

N2HDM exhibits a different vacuum structure compared to 2HDM (for a detailed
review of the vacuum structure of the N2HDM, please refer to [47]). General
properties of vacuum structure has been studied in [48]). As the loop
corrections and finite temperature effects may change the vacuum
configuration, the doublets and singlet can be expressed in the general vacuum
structure as:

$\Phi_{1}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\rho_{1}+i\eta_{1}\\\
\left(\zeta_{1}+\omega_{1}+i\psi_{1}\right)\end{array}\right),\quad\Phi_{2}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\rho_{2}+\omega_{CB}+i\eta_{2}\\\
\left(\zeta_{2}+\omega_{2}+i(\psi_{2}+\omega_{CP})\right)\end{array}\right),\quad\Phi_{S}=\zeta_{3}+\omega_{S},$
(3)

The Higgs fields are expanded in terms of the charged sectors $\rho_{i}$ and
$\eta_{i}$ ($i=1,2$) and the neutral CP-even and CP-odd fields $\zeta_{i}$ and
$\psi_{i}$ respectively. Our objective is to secure a minimum that conserves
both electric charge and CP while yielding three CP-even massive scalars. This
requirement, from the standpoint of vacuum structure suggests $\omega_{1}\neq
0,~{}\omega_{2}\neq 0~{}\rm{and}~{}\omega_{s}\neq 0$ while the CB and the CP
components are zero. This in turn, provides us (2). Because of the uniqueness
of the vacuum structure, the tree level potential itself contributes to the
strength of the electroweak phase transition by enhancing it. This, in turn,
increases the energy density of the stochastic GWs produced due to this
phenomenon as we shall see in the later part of this paper.

##  3 The One Loop Finite Temperature Effective potential

The effective potential $V_{\mathrm{eff}}$, encapsulates the influence of
radiative corrections on the scalar potential within the theory, particularly
at zero temperature [49]. At the one-loop order, $V_{\mathrm{eff}}$ can be
expressed as the sum of the N2HDM tree-level potential $V_{\mathrm{tree}}$ (2)
and the Coleman-Weinberg potential $V_{\rm CW}$. The Coleman-Weinberg
potential is well-established and is defined by:

$V_{\text{CW}}(\phi_{i})=\sum_{j}\frac{n_{j}}{64\pi^{2}}(-1)^{2s_{i}}\,m_{j}(\phi_{i})^{4}\left[\ln\left(\frac{|m_{j}(\phi_{i})^{2}|}{\mu^{2}}\right)-c_{j}\right].$
(4)

Here, $m_{j}(\phi_{i})$ represents the field-dependent tree-level mass of the
particle species $j$ in our model, where $n_{j}$ denotes its corresponding
number of degrees of freedom and $s_{j}$ represents the particle spin. To
compute the zero-temperature loop-corrected scalar masses, we incorporate
counter terms into $V_{\mathrm{eff}}$, yielding:

$V_{\rm CT}=\sum_{i}\frac{\partial V_{0}}{\partial p_{i}}\delta
p_{i}+\sum_{k}(\phi_{k}+v_{k})\delta T_{k}\ ,$ (5)

where $p_{i}$ represents the parameters of the tree-level potential.

###  3.1 Finite Temperature Effective Potential

To explore the thermal evolution of the N2HDM, it’s essential to compute the
effective potential, accounting for finite temperature corrections. The one-
loop effective potential, now written as $V_{\rm
tot}(\Phi_{1},\Phi_{2},\Phi_{S},T)$, under finite temperature corrections is
expressed as, [50]:

$V_{\rm tot}(\Phi_{1},\Phi_{2},\Phi_{S},T)\equiv V_{\rm tree}+V_{\rm
CW}+V_{T}\,,$ (6)

where $V_{T}$ is the one-loop thermal potential,

$V_{T}(\phi_{i})=\sum_{j}\,\frac{n_{j}\,T^{4}}{2\pi^{2}}\,J_{\pm}\left(\frac{m_{j}^{2}(\phi_{i})}{T^{2}}\right),$
(7)

with the thermal integrals for fermionic ($J_{+}$) and bosonic ($J_{-}$)
sectors:

$J_{\pm}\left(\frac{m_{j}^{2}(\phi_{i})}{T^{2}}\right)=\mp\int_{0}^{\infty}dx\,x^{2}\,\log\left[1\pm\exp\left(-\sqrt{x^{2}+\frac{m_{j}^{2}(\phi_{i})}{T^{2}}}\right)\right]\,,$
(8)

which vanish as $T\rightarrow 0$ (assuming $m^{2}_{j}$ is positive). In the
high temperature limit, it is convenient to expand the thermal functions
$J_{\pm}$:

$\displaystyle J_{-}(y)$ $\displaystyle\approx$
$\displaystyle-\frac{\pi^{4}}{45}+\frac{\pi^{2}}{12}y-\frac{\pi}{6}y^{\frac{3}{2}}-\frac{1}{32}\,y^{2}\,\text{log}\left(\frac{|y|}{a_{b}}\right)+\mathcal{O}(y^{3})\,,$
$\displaystyle J_{+}(y)$ $\displaystyle\approx$
$\displaystyle-\frac{7\pi^{4}}{360}+\frac{\pi^{2}}{24}y+\frac{1}{32}\,y^{2}\,\text{log}\left(\frac{|y|}{a_{f}}\right)+\mathcal{O}(y^{3})\,\quad\quad\mathrm{for}\,\left|y\right|\ll
1\,,$ (9)

with $a_{b}=\pi^{2}\text{exp}(3/2-2\gamma_{E})$ and
$a_{f}=16\pi^{2}\text{exp}(3/2-2\gamma_{E})$, $\gamma_{E}=0.57721\ldots$ being
the Euler-Mascheroni constant. In our analysis, we focus solely on the
contribution from the top quark, as the impact of other fermions within the SM
is relatively minor due to their small Yukawa couplings. Additionally, when
the field-dependent masses at a given temperature are significantly smaller
than the temperature of the plasma, i.e., $\frac{m_{i}^{2}}{T^{2}}<<1$, the
thermal function allows for a high-temperature expansion. This feature proves
highly beneficial for practical applications.

In another scenario, when considering a given plasma temperature, if the
field-dependent mass of a particle surpasses the temperature significantly,
i.e., $\frac{m_{i}(\Phi_{i})^{2}}{T^{2}}>>1$, both bosonic and fermionic
thermal functions exhibit behavior akin to exponentially decreasing functions
[50]. Consequently, in the thermal effective potential, the contribution of a
particle with a field-dependent mass greater than the temperature becomes
nearly negligible. The initial terms on the right-hand side of equations (9)
are devoid of classical background fields, rendering them inconsequential in
computing the critical temperature, $T_{c}$.

##  4 Stochastic Gravitational Waves

###  4.1 First Order Electroweak Phase Transition and Entropy Production

The study of electroweak phase transition (EWPT) is important in the study of
electroweak baryogenesis (EWBG) and GWs. A model which undergoes a strong
first order EWPT can successfully create EWBG and produce stochastic GW during
the PT. A strong FOPT is necessary to suppress processes such as $SU(2)_{L}$
sphalerons, which could potentially erase the baryon asymmetry after its
production. A strong first order EWPT can result in the production and release
of entropy into the primeval plasma, which dilutes any preexisting baryon
asymmetry and the energy density of the GWs. Hence, at this point it is
necessary to look into the entropy production due to EWPT in the framework of
N2HDM.

It is an established fact that the entropy density within the primordial
plasma remains constant throughout cosmological expansion, provided the plasma
remains in a state of thermal equilibrium with negligible chemical potentials
for all particle species. It can be expressed as:

$s=\frac{\mathcal{P}+\rho}{T}a^{3}=\rm{const,}$ (10)

where $\mathcal{P}~{}\rm{and}\rho$ are the pressure and the energy density and
$T$ is the temperature of the plasma while $a(t)$ is the cosmological scale
factor.

The Lagrangian density for the electroweak theory in N2HDM can be expressed
as:

$\mathcal{L}=\mathcal{L}_{f}+\mathcal{L}_{\rm Yuk}+\mathcal{L}_{\rm
gauge,kin}+\mathcal{L}_{\rm Scalar}.$ (11)

In (11), $\mathcal{L}_{f}$ is the kinetic term which originates from the
fermionic sector of the model and is given by:

$\displaystyle\mathcal{L}_{f}$ $\displaystyle=$
$\displaystyle\sum_{j}i(\bar{\Psi}^{(j)}_{L}\not{D}\Psi^{(j)}_{L}+\bar{\Psi}^{(j)}_{R}\not{D}\Psi^{(j)}_{R})$
(12) $\displaystyle=$ $\displaystyle
i\bar{\Psi}_{L}\gamma^{\xi}(\partial_{\mu}+igW_{\mu}+ig^{\prime}Y_{L}B_{\mu})\Psi_{L}$
$\displaystyle+i\bar{\Psi}_{R}\gamma^{\mu}(\partial_{\mu}+igW_{\mu}+ig^{\prime}Y_{R}B_{\mu})\Psi_{R}.$

Here $L$ and $R$ denote the left and right chiral fields of the fermion,
respectively. $\not{D}$ represents the gauge covariant derivative in Feynman
notation, and $j$ encompasses all fermionic species (i.e., the field
$\Psi_{j}$). The symbol $g$ represents the coupling constant. The partial
derivative goes as $\partial_{0}=d/dt$ and $\mu$ runs from $0$ to $3$, [51].
The second term on the R.H.S of (11) originates from the Yukawa interaction
and is given by:

$\displaystyle\mathcal{L}_{\rm
Yuk}=-\left[y_{e}\bar{e_{R}}\Phi_{a}^{\dagger}L_{L}+y_{e}^{*}\bar{L_{L}}\Phi_{a}^{\dagger}e_{R}+\cdots\right],$
(13)

where $y_{e}$ is dimensionless and complex constant. $\Phi_{a}$ (where
$a=1,2$) represents an $SU(2)_{L}$ doublet, and to maintain gauge invariance,
it’s coupled with another $SU(2)_{L}$ fermion, $L_{L}$. Similarly, $e_{R}$ and
other fermions such as quarks and neutrinos are right chiral fields,
consistent with the aforementioned coupling. $\mathcal{L}_{\rm gauge,\,kin}$
represents $U(1)$ invariant kinetic term of four gauge bosons
($W^{i},\,i=1,2,3$, and $B$). It can be written as:

$\displaystyle\mathcal{L}_{\rm
gauge,\,kin}=-\frac{1}{4}G^{i}_{\mu\nu}{G^{i}}^{\mu\nu}-\frac{1}{4}F^{B}_{\mu\nu}{F^{B}}^{\mu\nu},$
(14)

where
$G^{i}_{\mu\nu}=\partial_{\mu}W^{i}_{\nu}-\partial_{\nu}W^{i}_{\mu}-g\epsilon^{ijk}W_{\mu}^{j}W_{\nu}^{k}$
and $F^{B}_{\mu\nu}=\partial_{\mu}B_{\nu}-\partial_{\nu}B_{\mu}$ with
$\epsilon_{ijk}$ is the Levi-Civita symbol. The last term in (11) deals with
the scalar sector which includes the two Higgs doublet and the scalar singlet
particle. It can be written as:

$\displaystyle\mathcal{L}_{\rm Scalar}$ $\displaystyle=$
$\displaystyle(D^{\mu}\Phi_{1})^{\dagger}(D_{\mu}\Phi_{2})+(D^{\mu}\Phi_{1})^{\dagger}(D_{\mu}\Phi_{2})+2((D^{\mu}\Phi_{1})^{\dagger}(D_{\mu}\Phi_{S})+(D^{\mu}\Phi_{2})^{\dagger}(D_{\mu}\Phi_{S}))$
(15) $\displaystyle-V_{\rm tot}(\Phi_{1},\Phi_{2},\Phi_{S},T),$

where

$V_{\rm tot}(\Phi_{1},\Phi_{2},\Phi_{S},T)\equiv V_{\rm tree}+V_{\rm
CW}+V_{T},$ (16)

where all the individual components of the RHS are mentioned in (2), (4) and
(7).

As the temperature of the universe decreases, a secondary local minimum
emerges. When the plasma temperature reaches $T_{c}$, this secondary local
minimum adopts the form
$(\left<\Phi_{1}\right>=v_{1},\left<\Phi_{2}\right>=v_{2},\left<\Phi_{S}\right>=v_{S})$
and becomes degenerate with the global minimum at
$(\left<\Phi_{1}\right>=\left<\Phi_{2}\right>=\left<\Phi_{S}\right>=0)$. The
structure of the potential after the symmetry breaking is shown in figure 1.
The critical temperature is determined by the following expression:

$V_{\rm tot}\left(\Phi_{1}=0,\Phi_{2}=0,\Phi_{S}=0,T_{c}\right)=V_{\rm
tot}\left(\Phi_{1}=v_{1},\Phi_{2}=v_{2},,\Phi_{S}=v_{S},T_{c}\right).$ (17)

Figure 1: Strong EWPT can be seen in this case when the symmetry breaking
takes place. The plot is made with respect to $\Phi_{S}$ as it contains the
cubic term in the potential in (2). The parameter space of BP1 (1) has been
followed here.

Assuming the field to be homogeneous, the energy density can be written as:

$\displaystyle\rho=$
$\displaystyle\partial^{0}\Phi_{a}^{\dagger}\partial^{0}\Phi_{a}-(\mathcal{W}^{0}\Phi_{a})^{\dagger}\mathcal{W}^{0}\Phi_{a}-(\mathcal{W}^{j}\Phi_{a})^{\dagger}\mathcal{W}_{j}\Phi_{a}$
$\displaystyle+\left[V_{\rm
tot}(\Phi_{1},\Phi_{2},\Phi_{S},T)-\mathcal{L}_{\rm
gauge,kin}-\mathcal{L}_{f}-\mathcal{L}_{\rm Yuk},\right]$ (18)

where $\mathcal{W}_{\mu}=gT^{i}W^{i}_{\mu}+g^{\prime}YB_{\mu}$. Because of the
homogeneity and isotropy condition the spatial derivatives of the Higgs fields
vanishes. Following (10), the sum of pressure density and energy density of
the plasma is given by:

$\displaystyle\rho+\mathcal{P}=2\partial^{0}\Phi_{a}\partial^{0}\Phi_{a}^{\dagger}-i({\mathcal{W}^{0}}\Phi_{a})^{\dagger}(\partial_{0}\Phi_{a})+i(\partial^{0}{\Phi_{a}}^{\dagger})\mathcal{W}_{0}\Phi_{a}.$
(19)

Near the electroweak phase transition, $g_{*}\sim 110$ represents the
effective number of relativistic particles. This value varies with temperature
and diminishes as the cosmological cooling proceeds. (18) remains valid in the
scenario of instantaneous thermalization. Hence it can be rewritten as:

$\displaystyle\rho=\dot{\Phi}_{a,{\rm min}}^{2}+V_{\rm
tot}(\Phi_{1},\Phi_{2},\Phi_{S}T)+\frac{g_{*}\pi^{2}}{30}T^{4}.$ (20)

The final term on the right-hand side of (20) accounts for the energy density
of relativistic particles that have yet to acquire mass by the time of the
electroweak phase transition (EWPT). This includes contributions from Yukawa
interactions between fermions and Higgs bosons, as well as the energy density
associated with fermions, gauge bosons, and the interaction between Higgs
bosons and gauge bosons.

Before the onset of EWPT, $T\gg\left(\mu_{\rm chm.\,pot.}-m\right)$, where
$\mu_{\rm chm.\,pot.}$ is the chemical potential of the particle species and
it is assumed to be negligible for fermions. The entropy conservation law for
the relativistic particles is:

$\displaystyle
s\left(t\right)=\frac{2}{45}\pi^{2}g_{\star,s}\,T\left(t\right)^{3},$ (21)

where $g_{\star,s}$ is effective degrees of freedom in the entropy which
varies with time and temperature. Although it’s influenced by the constituents
of the primordial hot mixture, the effective number of relativistic degrees of
freedom, denoted as $g_{\star}$, might not always align with the count of
relativistic species. This discrepancy is discussed further in [52]. However,
for the purposes of this project, we assumed $g_{\star,s}\approx g_{\star}$.
From (10) and (21), we derive $T\sim a^{-1}$, indicating that temperature
scales inversely with the scale factor. As the universe expands and cools, it
reaches a state of thermal disequilibrium at some stage. Consequently, the
value of $s$ and, correspondingly, $g_{\star}\left(T\right)a^{3}T^{3}$ may
have increased, as entropy can either rise or remain constant.

As the temperature decreases from $T_{c}$ during EWPT, specific components of
the relativistic plasma reach their decoupling temperatures. At this point,
these components transition from being relativistic to non-relativistic and
acquire mass. The temperature at which decoupling occurs hinges upon the
masses of the components and their respective coupling constants. This
decoupling process leads to alterations in the effective number of
relativistic degrees of freedom, $g_{\star}$, within the plasma. EWPT presents
a potential scenario for Electroweak Baryogenesis to occur, establishing
thermal disequilibrium in the universe according to Sakharov’s principle [1].

As a consequence, the law of entropy conservation is violated, resulting in a
net increase in entropy. This increase in entropy is primarily driven by the
decoupling process, leading to a change in the effective number of
relativistic degrees of freedom, $g_{\star}$, within the relativistic plasma.
The predominant contribution to this increase in entropy stems from the
presence of the heaviest particle with a mass $m(T)$ below the prevailing
temperature $T$. Consequently, within this temperature range, the product
$g_{\star}(T)T$ remains constant, causing the entropy to rise proportionally
to $a^{3}T^{3}$. Since the final temperature $T=m(T_{f})$ (below which a new
particle begins to dominate) is independent of $g_{\star}$, we observe an
increase in entropy against the backdrop of $g_{\star}a^{3}T^{3}$.

For brevity, we assume the oscillations of $\phi$ around $\phi_{min}$ undergo
rapid damping. Consequently, we set $\dot{\phi}=\dot{\phi}_{min}$ and
disregard higher-order terms, as the evolution of $\phi_{min}$ is primarily
driven by the gradual expansion of the universe, which occurs at a much slower
pace. This simplification results in a single differential equation for the
temperature (or scale factor). In principle, this can be numerically computed
by solving the relevant Klein-Gordon equation, accounting for damping effects
caused by particle production.

To compute the entropy production resulting from the EWPT, which subsequently
impacts the amplitude of stochastic gravitational waves generated during the
EWPT, it’s essential to solve the evolution equation governing the
conservation of energy density,

$\dot{\rho}=-3\mathcal{H}(\rho+\mathcal{P}).$ (22)

In the subsequent analysis, we utilize the cosmoTransitions package, [53],
which specializes in handling a range of properties and functionalities
associated with phase transitions. This package is employed to compute
essential parameters such as the critical temperature ($T_{c}$) and the vacuum
expectation value ($v_{c}$), as well as $V_{\rm eff}(T)$, for each of the
benchmark points under consideration. In order to achieve a first order phase
transition the ratio $v_{c}/T_{c}\gg 1$. The set of benchmark points (BPs) are
shown in Table 1

BPs |  $\lambda_{1}$ |  $\lambda_{2}$ |  $\lambda_{3}$ |  $\lambda_{4}$ |  $\lambda_{5}$ |  $a_{1}$ |  $a_{2}$ |  $a_{3}$ |  $\lambda_{6}$ |  $\lambda_{7}$ |  $a_{4}$ |  $\lambda_{8}$
---|---|---|---|---|---|---|---|---|---|---|---|---
BP1 |  1.28 |  1.28 |  0.7 |  0.37 |  0.37 |  0.5 |  0.5 |  0.1 |  0.5 |  0.5 |  0.1 |  0.5
BP2 |  1.28 |  0.002 |  0.7 |  0.24 |  0.24 |  0.7 |  0.7 |  0.1 |  1.5 |  1.0 |  0.1 |  0.5
BP3 |  0.263 |  0.258 |  1.06 |  0.3 |  1.35 |  0.1 |  0.5 |  0.5 |  1.35 |  1.35 |  0.5 |  1.35
BP4 |  1.14 |  0.37 |  0.5 |  0 |  0 |  1.2 |  1.2 |  0.8 |  1.4 |  1.4 |  0.8 |  0.65
BP5 |  4.13 |  0.22 |  4.15 |  0 |  0 |  1.8 |  1.8 |  1.2 |  1.6 |  1.6 |  1.2 |  1.6
Table 1: Estimates of the various parameters in (2). Please note that the
$\lambda_{i},~{}(i=1-8))$ are dimensionless constants while
$a_{i},~{}(i=1-4))$ are in GeV.

After numerically solving (22) for the various BPs in Table 1, we have plotted
out both the evolution of entropy as a function of the scale factor $a/a_{c}$
and the entropy release (or the net increase in entropy production) for each
BPs have been calculated and plotted in figure 2. The values of $T_{c}$ and
the entropy release is provided in Table 2.

Figure 2: The figures show the evolution of entropy, which is made dimensionless for brevity as a function of the ratio of scale factor to the scale factor corresponding to the transition temperature for all different BPs presented in Table 1. BPs |  $T_{c}(\rm{GeV})$ |  $\frac{\delta s}{s}($%$)$
---|---|---
BP1 |  129 |  30
BP2 |  96 |  40
BP3 |  119 |  34
BP4 |  125 |  33
BP5 |  89 |  44
Table 2: Values of the critical/transition temperature $T_{c}$ and entropy
production corresponding to BPs in Table 1.

We can see from Figure 2 that during the time of the EWPT, there is a
significant increase in the entropy of the plasma. The entropy starts
increasing from the onset of the PT and when the PT comes to an end, it
becomes a constant value. By constant value, it is meant that the entropy
curve follows the regular $~{}a^{-3}$ nature, with no change in the slope of
the curve. In terms of the net entropy released, it can be said that during
the phase transition, the entropy density increased by a certain amount,
provided in Table 2 and when the period of phase transition ended, it became
constant. It is evident from Table 2 that stronger the phase transition, the
larger is the production of entropy. In the Table 2, it can be seen that the
strength of the phase transition is strongest for BP5, which corresponds to
$\sim 44\%$ release in entropy compared to BP1, which has the least strength
among the $5$ BPs, which corresponds to $\sim 30\%$ increase in the entropy of
the plasa after the EWPT. The primary factor behind this surplus in entropy
production stems from the additional scalar sector generated in the N2HDM,
which notably contribute to this phenomenon. Notably, contributions from
lighter particles such as electrons and neutrinos mirror those observed in the
SM. In the next section we will show how the strength of the phase transition
is related to the production of stochastic GW.

###  4.2 Stochastic Gravitational Waves

A first-order phase transition occurring in the early universe serves as a
potential source of stochastic GW. Typically, this transition happens roughly
$\sim 10^{-11}~{}\rm{sec}$ after the big bang, preceding the initiation of Big
Bang Nucleosynthesis. As the temperature of the universe decreases, it
undergoes a first-order phase transition at $T=T_{c}$. Above this critical
temperature, the universe’s symmetry is intact. However, as the temperature
drops below $T_{c}$, a second degenerate minimum emerges, signaling the onset
of a first-order phase transition. During this transition, latent heat is
released, with a portion contributing to the generation of GWs, while the
remainder is absorbed by the plasma. Bubble nucleation governs the progression
of a FOPT, making the calculation of the nucleation temperature pivotal for
determining essential phenomenological parameters crucial for estimating
gravitational wave spectra. As nucleation occurs at temperatures below
$T_{c}$, the probability of tunneling $\Gamma(T)$ from the false vacuum to the
true one can be expressed as given in [54],

$\Gamma(T)\approx T^{4}\left(\frac{S_{E}}{2\pi
T}\right)^{3/2}\;e^{-\frac{S_{E}}{T}},$ (23)

where $S_{E}$ is the Euclidean action, otherwise known as the bounce action.
It is given by, [55]:

$\displaystyle S_{E}=\int_{0}^{\infty}4\pi
r^{2}dr\left(V_{T}(\phi,T)+\frac{1}{2}\;\left(\frac{d\phi(r)}{dr}\right)^{2}\right).$
(24)

Here $\phi$ is the scalar dynamical field whose classical solution can be
represented in terms of the radial coordinate $r$ as:

$\frac{d^{2}\phi}{dr^{2}}+\frac{2}{r}\frac{d\phi}{dr}=\frac{dV_{T}(\phi,T)}{dr},$
(25)

More details can be found in [56, 57]. The necessary boundary conditions to
solve (25) are:

$\frac{d\phi}{dr}=0~{}\rm{when}~{}r\rightarrow 0,$ (26)

and

$\phi(r)\rightarrow\phi_{\rm false}~{}\rm{when}~{}r\rightarrow\infty.$ (27)

In (27) $\phi_{\rm false}$ denotes the four dimensional field values at the
false vacuum. The package cosmoTransitions has been used to solve (25). The
key parameters crucial for estimating gravitational wave spectra originating
from first-order phase transitions are the relative change in energy density
during the phase transition ($\alpha$) and the inverse of the duration of the
phase transition ($\beta$). Both $\alpha$ and $\beta$ are defined at the
nucleation temperature $T_{n}$. The parameter $\alpha$ is calculated as:

$\alpha=\frac{\Delta\rho}{\rho_{\rm rad}},$ (28)

where $\rho_{\rm rad}$ is the radiation energy density and $\Delta\rho$ is the
released latent heat released during the phase transition, [58]. It is given
by, [59]:

${\Delta\rho=\left[V_{T}(\phi_{0},T)-T\frac{dV_{T}(\phi_{0},T)}{dT}\right]_{T=T_{n}}-\left[V_{T}(\phi_{n},T)-T\frac{dV_{T}(\phi_{n},T)}{dT}\right]_{T=T_{n}}.}$
(29)

$\phi_{0}$ and $\phi_{n}$ are the field values at false and true vacuum and
$V_{T}(\phi,T)$ is the finite temperature effective potential and in our case
it is same as (6). The $\beta$ parameter is defined as:

$\frac{\beta}{H_{\ast}}=T\frac{d}{dT}\left(\frac{S_{E}}{T}\right)\Bigg{|}_{T=T_{\ast}}\equiv\,T\frac{d}{dT}\left(\frac{S_{E}}{T}\right)\Bigg{|}_{T=T_{n}},$
(30)

where $H_{*}$ is the expansion rate of the Universe during the PT and $T_{*}$
stands for the PT temperature. In our analysis, we have assumed $T_{*}=T_{n}$,
where $T_{n}$ is the nucleation temperature. The various values of $\alpha$
and $\beta/H_{*}$ are shown in Table 3 for various benchmark points.

BPs |  $\alpha$ |  $\beta/H_{*}$
---|---|---
BP1 |  0.05 |  8930.8
BP2 |  0.4 |  1202.5
BP3 |  0.6 |  13866.3
BP4 |  0.1 |  1295.2
BP5 |  0.75 |  1757.2
Table 3: Estimates of the parameters $\alpha$ and $\beta$

There are three primany components which contributes to the energy density of
the stochastic GW. They are the contributions from the bubble wall collisions,
the sound waves and the magneto-hydrodynamic turbulence (MHD). Thus the net
energy spectrum of stochastic GW can be approximated as the sum of all the
three components:

$\Omega_{\rm GW}h^{2}\approx\Omega_{\rm col}h^{2}+\Omega_{\rm
sw}h^{2}+\Omega_{\rm tur}h^{2},$ (31)

where the first term on the RHS of (31) is the contribution from the bubble
wall collisions, the second term corresponds to the sound-waves and the third
term gives the contribution from MHD, [60]. Also,
$h=H_{0}/(100\text{\,km}\cdot\text{sec}^{-1}\cdot\text{Mpc}^{-1})$, where
$H_{0}$ is the value of the Hubble constant at the present day universe, see
[61]. The portion of the overall gravitational wave energy density stemming
from bubble wall collisions can be determined utilizing the envelope
approximation, [62], and can be expressed a function of frequency $f$ as:

$\Omega_{\rm col}h^{2}=1.67\times
10^{-5}{\left(\frac{\beta}{H_{\ast}}\right)^{-2}}\left(\frac{\kappa_{c}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g^{\ast}}\right)^{1/3}\left(\frac{0.11v^{3}_{w}}{0.42+v^{2}_{w}}\right)\frac{3.8\left(f/f_{\rm
col}\right)^{2.8}}{1+2.8\left(f/f_{\rm col}\right)^{3.8}}\;.$ (32)

Here $v_{w}$ is the bubble wall velocity and $\kappa_{c}$ is the efficiency
factor of bubble collision. It is given by:

$\kappa_{c}=\frac{0.715\alpha+\frac{4}{27}\sqrt{\frac{3\alpha}{2}}}{1+0.715\alpha}.$
(33)

$f_{col}$ is the red-shifted frequency peak and can be written as:

$f_{\rm col}=16.5\times
10^{-6}\left(\frac{f_{\ast}}{\beta}\right)\left(\frac{\beta}{H_{\ast}}\right)\left(\frac{T_{n}}{100{\rm\,GeV}}\right)\left(\frac{g^{\ast}}{100}\right)^{1/6}{\rm
Hz},$ (34)

where $f_{\ast}/\beta$ is called the fitting function and it is given by:

$\frac{f_{\ast}}{\beta}=\frac{0.62}{1.8-0.1v_{w}+v^{2}_{w}}.$ (35)

In what follows we will assume $v_{w}=1$. This is because expanding bubbles
can obtain relativistic terminal velocity. The contribution to the energy
density of GW from the sound waves is given by, [63, 64, 65]:

$\Omega_{\rm sw}h^{2}=2.65\times 10^{-6}\;\Upsilon(\tau_{\rm
sw})\left(\frac{\beta}{H_{\ast}}\right)^{-1}v_{w}\left(\frac{\kappa_{\rm
sw}\alpha}{1+\alpha}\right)^{2}\left(\frac{g^{\ast}}{100}\right)^{1/3}\left(\frac{f}{f_{\rm
sw}}\right)^{3}\left[\frac{7}{4+3\left(f/f_{\rm sw}\right)^{2}}\right]^{7/2},$
(36)

where $\kappa_{\rm sw}$ denotes the efficiency factor for the sound wave
contribution, indicating the latent heat transformed into bulk motion of the
plasma, thereby emitting gravitational waves. This expression holds true in
the limit as the wall velocity approaches $v_{w}\to 1$, and can be expressed
as:

$\kappa_{\rm
sw}\simeq\left[\frac{\alpha}{0.73+0.083\sqrt{\alpha}+\alpha}\right].$ (37)

The peak frequency for the contribution of sound waves is denoted by $f_{\rm
sw}$ and is given by:

$f_{\rm sw}=1.9\times
10^{-5}\left(\frac{1}{v_{w}}\right)\left(\frac{\beta}{H_{\ast}}\right)\left(\frac{T_{n}}{100{\rm~{}GeV}}\right)\left(\frac{g^{\ast}}{100}\right)^{1/6}{\rm
Hz}.$ (38)

The suppression of the contribution of the sound waves to the total energy
density of the GW energy spectrum is caused by $\Upsilon(\tau_{\rm sw})$,
which comes into play due to the finite lifetime of the sound waves and it is
given by:

$\Upsilon(\tau_{\rm sw})=1-\frac{1}{\sqrt{1+2\tau_{\rm sw}H_{\ast}}},$ (39)

with $\tau_{\rm sw}$ being the lifetime of the sound waves. The contribution
from the MHD, which arises due to the complete ionization of the plasma, [66],
is given by:

$\Omega_{\rm tur}h^{2}=3.35\times
10^{-4}\left(\frac{\beta}{H_{\ast}}\right)^{-1}v_{w}\left(\frac{\kappa_{\rm
tur}\alpha}{1+\alpha}\right)^{3/2}\left(\frac{100}{g^{\ast}}\right)^{1/3}\left[\frac{\left(f/f_{\rm
tur}\right)^{3}}{\left[1+\left(f/f_{\rm
tur}\right)\right]^{11/3}\left(1+\frac{8\pi f}{h_{\ast}}\right)}\right],$ (40)

where
$h_{\ast}=16.5\times\left(\frac{T_{n}}{100{\rm~{}GeV}}\right)\left(\frac{g^{\ast}}{100}\right)^{1/6}{\rm
Hz}$, the inverse Hubble time during GW production, red-shifted to today. The
peak frequency $f_{\rm tur}$ is given by,

$f_{\rm tur}=2.7\times
10^{-5}\frac{1}{v_{w}}\left(\frac{\beta}{H_{\ast}}\right)\left(\frac{T_{n}}{100{\rm~{}GeV}}\right)\left(\frac{g^{\ast}}{100}\right)^{1/6}{\rm
Hz.}$ (41)

$\kappa_{\rm tur}$ is defined as $\kappa_{\rm tur}=\epsilon\kappa_{\rm{sw}}$,
where $\epsilon$ stands for the fraction of the bulk motion which is turbulent
and previous results, [66], suggests $\kappa_{\rm tur}=0.1\kappa_{\rm sw}$,
which is used in the calculations.

Figure 3 shows the plots of the total energy density of the produced
stochastic GWs (equivalent to (31)), considering all the aforementioned
points. The produced results has been compared with the predictions from
various GW interferometers LISA, BBO, DECIGO and others [67, 68, 69, 70, 71,
72, 73, 74, 75, 76, 77, 78].

Figure 3: GW spectrum for various benchmark points corresponding to Table 3.
The results of our calculations are compared with the available sensitivity
plots provided by the various GW interferometers for future observations. [67,
68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78]

The parameter $\alpha$ is directly related to the energy released during the
phase transition. Therefore, a stronger phase transition should result in a
higher value of $\alpha$. It is obvious from the plot in Figure 3 that the
strength of the phase transition corresponding to BP5 is the strongest and
hence we see that the amplitude of the GW spectrum has the highest value when
BP5 is considered. Similarly, for BP1, we find the strength of the phase
transition is the weakest because $\alpha$ has the lowest value corresponding
to BP1 and hence this is reflected on the amplitude of the GW spectrum.

##  5 Conclusion

In this work, we studied two implications of a strong first order EWPT within
the framework of the next-to-minimal Two Higgs Doublet Model. We considered
only one $\mathbb{Z}_{2}$ symmetry imposed on the Higgs doublets and the
singlet, thus allowing the singlet to have cubic terms in the potential which
gives rise to a strong FOPT inherently, as can be seen in Figure 1. We
proceeded to calculate the entropy production due to this EWPT and we found
out that there is a significant increase/influx of entropy into plasma during
the phase transition. Furthermore, we have studied the production of
stochastic gravitational waves which can be produced during the EWPT and we
indeed found prominent signatures of GW signals and they overlap well with the
sensitivity curves provided by the GW interferometers and are open to
detection in the near future.

We saw how the strength of the phase transition effects the entropy production
and the production of GW. A stronger PT results in a higher production of both
entropy and GW. However, in our calculations, we have assumed the bubble wall
velocity $v_{w}\sim 1$. But a precise calculation and determination of this
bubble wall velocity, even though, non-trivial, is possible as shown in [79,
80, 81]. This might have an effect on the net production of GW but analysis of
(31) considering non-unitary values of $v_{w}$ is beyond the scope of this
paper.

##  6 Acknowledgement

The work of AC was supported by the Japan Society for the Promotion of Science
(JSPS) as a part of the JSPS Postdoctoral Program (Standard), grant number:
23KF0289. The work of KK was supported by the KAKENHI grant No. 24H01825. AC
also thanks Kousik Loho for useful discussion.

## References

  * [1] A. D. Sakharov, Pisma Zh. Eksp. Teor. Fiz. 5, 32-35 (1967) doi:10.1070/PU1991v034n05ABEH002497
  * [2] S. Chatrchyan et al. [CMS], Phys. Lett. B 716, 30-61 (2012) doi:10.1016/j.physletb.2012.08.021 [arXiv:1207.7235 [hep-ex]].
  * [3] A. Chaudhuri and A. Dolgov, JCAP 01 (2018), 032 doi:10.1088/1475-7516/2018/01/032 [arXiv:1711.01801 [hep-ph]].
  * [4] K. Kajantie, M. Laine, K. Rummukainen and M. E. Shaposhnikov, Phys. Rev. Lett. 77, 2887-2890 (1996) doi:10.1103/PhysRevLett.77.2887 [arXiv:hep-ph/9605288 [hep-ph]].
  * [5] A. Chaudhuri and J. Das, Phys. Rev. D 106, no.9, 095016 (2022) doi:10.1103/PhysRevD.106.095016 [arXiv:2206.08699 [hep-ph]].
  * [6] D. E. Morrissey and M. J. Ramsey-Musolf, New J. Phys. 14 (2012), 125003 doi:10.1088/1367-2630/14/12/125003 [arXiv:1206.2942 [hep-ph]].
  * [7] M. Carena, Z. Liu and Y. Wang, JHEP 08, 107 (2020) doi:10.1007/JHEP08(2020)107 [arXiv:1911.10206 [hep-ph]].
  * [8] C. W. Chiang and B. Q. Lu, JHEP 07, 082 (2020) doi:10.1007/JHEP07(2020)082 [arXiv:1912.12634 [hep-ph]].
  * [9] Z. Kang, P. Ko and T. Matsui, JHEP 02, 115 (2018) doi:10.1007/JHEP02(2018)115 [arXiv:1706.09721 [hep-ph]].
  * [10] A. Mazumdar and G. White, Rept. Prog. Phys. 82, no.7, 076901 (2019) doi:10.1088/1361-6633/ab1f55 [arXiv:1811.01948 [hep-ph]].
  * [11] R. Apreda, M. Maggiore, A. Nicolis and A. Riotto, Nucl. Phys. B 631, 342-368 (2002) doi:10.1016/S0550-3213(02)00264-X [arXiv:gr-qc/0107033 [gr-qc]].
  * [12] C. Grojean, G. Servant and J. D. Wells, Phys. Rev. D 71, 036001 (2005) doi:10.1103/PhysRevD.71.036001 [arXiv:hep-ph/0407019 [hep-ph]].
  * [13] D. J. Weir, Phil. Trans. Roy. Soc. Lond. A 376, no.2114, 20170126 (2018) [erratum: Phil. Trans. Roy. Soc. Lond. A 381, 20230212 (2023)] doi:10.1098/rsta.2017.0126 [arXiv:1705.01783 [hep-ph]].
  * [14] J. Ellis, M. Lewicki and J. M. No, JCAP 04, 003 (2019) doi:10.1088/1475-7516/2019/04/003 [arXiv:1809.08242 [hep-ph]].
  * [15] T. Alanne, T. Hugle, M. Platscher and K. Schmitz, JHEP 03, 004 (2020) doi:10.1007/JHEP03(2020)004 [arXiv:1909.11356 [hep-ph]].
  * [16] C. Caprini, M. Chala, G. C. Dorsch, M. Hindmarsh, S. J. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No and K. Rummukainen, et al. JCAP 03, 024 (2020) doi:10.1088/1475-7516/2020/03/024 [arXiv:1910.13125 [astro-ph.CO]].
  * [17] A. Chaudhuri, M. Y. Khlopov and S. Porey, Symmetry 14, no.2, 271 (2022) doi:10.3390/sym14020271 [arXiv:2110.14161 [hep-ph]].
  * [18] Y. Kikuta, K. Kohri and E. So, [arXiv:1405.4166 [hep-ph]].
  * [19] F. Costa, S. Khan and J. Kim, JHEP 06, 026 (2022) doi:10.1007/JHEP06(2022)026 [arXiv:2202.13126 [hep-ph]].
  * [20] A. Ahriche, K. Hashino, S. Kanemura and S. Nasri, Phys. Lett. B 789, 119-126 (2019) doi:10.1016/j.physletb.2018.12.013 [arXiv:1809.09883 [hep-ph]].
  * [21] V. R. Shajiee and A. Tofighi, Eur. Phys. J. C 79, no.4, 360 (2019) doi:10.1140/epjc/s10052-019-6881-6 [arXiv:1811.09807 [hep-ph]].
  * [22] A. Alves, T. Ghosh, H. K. Guo, K. Sinha and D. Vagie, JHEP 04, 052 (2019) doi:10.1007/JHEP04(2019)052 [arXiv:1812.09333 [hep-ph]].
  * [23] A. Chatterjee, A. Datta and S. Roy, JHEP 06, 108 (2022) doi:10.1007/JHEP06(2022)108 [arXiv:2202.12476 [hep-ph]].
  * [24] W. Wang, K. P. Xie, W. L. Xu and J. M. Yang, Eur. Phys. J. C 82, no.12, 1120 (2022) doi:10.1140/epjc/s10052-022-11077-3 [arXiv:2204.01928 [hep-ph]].
  * [25] M. Kakizaki, S. Kanemura and T. Matsui, Phys. Rev. D 92, no.11, 115007 (2015) doi:10.1103/PhysRevD.92.115007 [arXiv:1509.08394 [hep-ph]].
  * [26] K. Hashino, M. Kakizaki, S. Kanemura and T. Matsui, Phys. Rev. D 94, no.1, 015005 (2016) doi:10.1103/PhysRevD.94.015005 [arXiv:1604.02069 [hep-ph]].
  * [27] K. Hashino, M. Kakizaki, S. Kanemura, P. Ko and T. Matsui, Phys. Lett. B 766, 49-54 (2017) doi:10.1016/j.physletb.2016.12.052 [arXiv:1609.00297 [hep-ph]].
  * [28] K. Hashino, M. Kakizaki, S. Kanemura, P. Ko and T. Matsui, JHEP 06, 088 (2018) doi:10.1007/JHEP06(2018)088 [arXiv:1802.02947 [hep-ph]].
  * [29] M. J. Ramsey-Musolf, JHEP 09, 179 (2020) doi:10.1007/JHEP09(2020)179 [arXiv:1912.07189 [hep-ph]].
  * [30] L. S. Friedrich, M. J. Ramsey-Musolf, T. V. I. Tenkanen and V. Tran, [arXiv:2203.05889 [hep-ph]].
  * [31] K. Hashino, R. Jinno, M. Kakizaki, S. Kanemura, T. Takahashi and M. Takimoto, Phys. Rev. D 99, no.7, 075011 (2019) doi:10.1103/PhysRevD.99.075011 [arXiv:1809.04994 [hep-ph]].
  * [32] G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher and J. P. Silva, Phys. Rept. 516, 1-102 (2012) doi:10.1016/j.physrep.2012.02.002 [arXiv:1106.0034 [hep-ph]].
  * [33] B. Grzadkowski and P. Osland, Phys. Rev. D 82, 125026 (2010) doi:10.1103/PhysRevD.82.125026 [arXiv:0910.4068 [hep-ph]].
  * [34] C. Y. Chen, M. Freid and M. Sher, Phys. Rev. D 89, no.7, 075009 (2014) doi:10.1103/PhysRevD.89.075009 [arXiv:1312.3949 [hep-ph]].
  * [35] A. Drozd, B. Grzadkowski, J. F. Gunion and Y. Jiang, JHEP 11, 105 (2014) doi:10.1007/JHEP11(2014)105 [arXiv:1408.2106 [hep-ph]].
  * [36] M. Muhlleitner, M. O. P. Sampaio, R. Santos and J. Wittbrodt, JHEP 03, 094 (2017) doi:10.1007/JHEP03(2017)094 [arXiv:1612.01309 [hep-ph]].
  * [37] G. C. Dorsch, S. J. Huber and J. M. No, JHEP 10, 029 (2013) doi:10.1007/JHEP10(2013)029 [arXiv:1305.6610 [hep-ph]].
  * [38] G. C. Dorsch, S. J. Huber, K. Mimasu and J. M. No, Phys. Rev. Lett. 113, no.21, 211802 (2014) doi:10.1103/PhysRevLett.113.211802 [arXiv:1405.5537 [hep-ph]].
  * [39] P. Basler, M. Krause, M. Muhlleitner, J. Wittbrodt and A. Wlotzka, JHEP 02, 121 (2017) doi:10.1007/JHEP02(2017)121 [arXiv:1612.04086 [hep-ph]].
  * [40] G. C. Dorsch, S. J. Huber, K. Mimasu and J. M. No, JHEP 12, 086 (2017) doi:10.1007/JHEP12(2017)086 [arXiv:1705.09186 [hep-ph]].
  * [41] J. Bernon, L. Bian and Y. Jiang, JHEP 05, 151 (2018) doi:10.1007/JHEP05(2018)151 [arXiv:1712.08430 [hep-ph]].
  * [42] S. Heinemeyer, C. Li, F. Lika, G. Moortgat-Pick and S. Paasch, Phys. Rev. D 106, no.7, 075003 (2022) doi:10.1103/PhysRevD.106.075003 [arXiv:2112.11958 [hep-ph]].
  * [43] T. Biekötter, M. Chakraborti and S. Heinemeyer, Eur. Phys. J. C 80, no.1, 2 (2020) doi:10.1140/epjc/s10052-019-7561-2 [arXiv:1903.11661 [hep-ph]].
  * [44] T. Biekötter, S. Heinemeyer, J. M. No, M. O. Olea and G. Weiglein, JCAP 06, 018 (2021) doi:10.1088/1475-7516/2021/06/018 [arXiv:2103.12707 [hep-ph]].
  * [45] I. P. Ivanov, Phys. Rev. D 77, 015017 (2008) doi:10.1103/PhysRevD.77.015017 [arXiv:0710.3490 [hep-ph]].
  * [46] A. Barroso, P. M. Ferreira and R. Santos, Phys. Lett. B 652, 181-193 (2007) doi:10.1016/j.physletb.2007.07.010 [arXiv:hep-ph/0702098 [hep-ph]].
  * [47] P. Basler, M. Mühlleitner and J. Müller, JHEP 05, 016 (2020) doi:10.1007/JHEP05(2020)016 [arXiv:1912.10477 [hep-ph]].
  * [48] P. M. Ferreira, M. Mühlleitner, R. Santos, G. Weiglein and J. Wittbrodt, JHEP 09, 006 (2019) doi:10.1007/JHEP09(2019)006 [arXiv:1905.10234 [hep-ph]].
  * [49] M. Quiros, Helv. Phys. Acta 67, 451-583 (1994)
  * [50] C. Delaunay, C. Grojean and J. D. Wells, JHEP 04, 029 (2008) doi:10.1088/1126-6708/2008/04/029 [arXiv:0711.2511 [hep-ph]].
  * [51] S. Karmakar and S. Rakshit, Phys. Rev. D 100, no.5, 055016 (2019) doi:10.1103/PhysRevD.100.055016 [arXiv:1901.11361 [hep-ph]].
  * [52] C. Lunardini and Y. F. Perez-Gonzalez, JCAP 08, 014 (2020) doi:10.1088/1475-7516/2020/08/014 [arXiv:1910.07864 [hep-ph]].
  * [53] C. L. Wainwright, Comput. Phys. Commun. 183, 2006-2013 (2012) doi:10.1016/j.cpc.2012.04.004 [arXiv:1109.4189 [hep-ph]].
  * [54] C. Grojean and G. Servant, Phys. Rev. D 75, 043507 (2007) doi:10.1103/PhysRevD.75.043507 [arXiv:hep-ph/0607107 [hep-ph]].
  * [55] A. D. Linde, Nucl. Phys. B 216, 421 (1983) [erratum: Nucl. Phys. B 223, 544 (1983)] doi:10.1016/0550-3213(83)90072-X
  * [56] I. Affleck, Phys. Rev. Lett. 46, 388 (1981) doi:10.1103/PhysRevLett.46.388
  * [57] A. D. Linde, Phys. Lett. B 70, 306-308 (1977) doi:10.1016/0370-2693(77)90664-5
  * [58] M. Kamionkowski, A. Kosowsky and M. S. Turner, Phys. Rev. D 49, 2837-2851 (1994) doi:10.1103/PhysRevD.49.2837 [arXiv:astro-ph/9310044 [astro-ph]].
  * [59] J. Kehayias and S. Profumo, JCAP 03, 003 (2010) doi:10.1088/1475-7516/2010/03/003 [arXiv:0911.0687 [hep-ph]].
  * [60] C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, A. Petiteau, P. Schwaller and G. Servant, et al. JCAP 04, 001 (2016) doi:10.1088/1475-7516/2016/04/001 [arXiv:1512.06239 [astro-ph.CO]].
  * [61] T. M. C. Abbott et al. [DES], Mon. Not. Roy. Astron. Soc. 480, no.3, 3879-3888 (2018) doi:10.1093/mnras/sty1939 [arXiv:1711.00403 [astro-ph.CO]].
  * [62] R. Jinno and M. Takimoto, Phys. Rev. D 95, no.2, 024009 (2017) doi:10.1103/PhysRevD.95.024009 [arXiv:1605.01403 [astro-ph.CO]].
  * [63] M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, Phys. Rev. Lett. 112, 041301 (2014) doi:10.1103/PhysRevLett.112.041301 [arXiv:1304.2433 [hep-ph]].
  * [64] M. Hindmarsh, Phys. Rev. Lett. 120, no.7, 071301 (2018) doi:10.1103/PhysRevLett.120.071301 [arXiv:1608.04735 [astro-ph.CO]].
  * [65] M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, Phys. Rev. D 96, no.10, 103520 (2017) [erratum: Phys. Rev. D 101, no.8, 089902 (2020)] doi:10.1103/PhysRevD.96.103520 [arXiv:1704.05871 [astro-ph.CO]].
  * [66] C. Caprini, R. Durrer and G. Servant, JCAP 12, 024 (2009) doi:10.1088/1475-7516/2009/12/024 [arXiv:0909.0622 [astro-ph.CO]].
  * [67] H. Kudoh, A. Taruya, T. Hiramatsu and Y. Himemoto, Phys. Rev. D 73, 064006 (2006) doi:10.1103/PhysRevD.73.064006 [arXiv:gr-qc/0511145 [gr-qc]].
  * [68] K. Yagi and N. Seto, Phys. Rev. D 83, 044011 (2011) [erratum: Phys. Rev. D 95, no.10, 109901 (2017)] doi:10.1103/PhysRevD.83.044011 [arXiv:1101.3940 [astro-ph.CO]].
  * [69] K. Nakayama and J. Yokoyama, JCAP 01, 010 (2010) doi:10.1088/1475-7516/2010/01/010 [arXiv:0910.0715 [astro-ph.CO]].
  * [70] A. Afzal et al. [NANOGrav], Astrophys. J. Lett. 951, no.1, L11 (2023) doi:10.3847/2041-8213/acdc91 [arXiv:2306.16219 [astro-ph.HE]].
  * [71] P. Amaro-Seoane et al. [LISA], [arXiv:1702.00786 [astro-ph.IM]].
  * [72] J. Aasi et al. [LIGO Scientific], Class. Quant. Grav. 32, 074001 (2015) doi:10.1088/0264-9381/32/7/074001 [arXiv:1411.4547 [gr-qc]].
  * [73] F. Acernese et al. [VIRGO], Class. Quant. Grav. 32, no.2, 024001 (2015) doi:10.1088/0264-9381/32/2/024001 [arXiv:1408.3978 [gr-qc]].
  * [74] T. Akutsu et al. [KAGRA], Nature Astron. 3, no.1, 35-40 (2019) doi:10.1038/s41550-018-0658-y [arXiv:1811.08079 [gr-qc]].
  * [75] J. Crowder and N. J. Cornish, Phys. Rev. D 72, 083005 (2005) doi:10.1103/PhysRevD.72.083005 [arXiv:gr-qc/0506015 [gr-qc]].
  * [76] S. Kawamura, M. Ando, N. Seto, S. Sato, T. Nakamura, K. Tsubono, N. Kanda, T. Tanaka, J. Yokoyama and I. Funaki, et al. Class. Quant. Grav. 28, 094011 (2011) doi:10.1088/0264-9381/28/9/094011
  * [77] M. Punturo, M. Abernathy, F. Acernese, B. Allen, N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsuglia and M. Beker, et al. Class. Quant. Grav. 27, 194002 (2010) doi:10.1088/0264-9381/27/19/194002
  * [78] D. Reitze, R. X. Adhikari, S. Ballmer, B. Barish, L. Barsotti, G. Billingsley, D. A. Brown, Y. Chen, D. Coyne and R. Eisenstein, et al. Bull. Am. Astron. Soc. 51, no.7, 035 (2019) [arXiv:1907.04833 [astro-ph.IM]].
  * [79] M. Dine, R. G. Leigh, P. Y. Huet, A. D. Linde and D. A. Linde, Phys. Rev. D 46, 550-571 (1992) doi:10.1103/PhysRevD.46.550 [arXiv:hep-ph/9203203 [hep-ph]].
  * [80] J. Ignatius, K. Kajantie, H. Kurki-Suonio and M. Laine, Phys. Rev. D 49, 3854-3868 (1994) doi:10.1103/PhysRevD.49.3854 [arXiv:astro-ph/9309059 [astro-ph]].
  * [81] G. D. Moore and T. Prokopec, Phys. Rev. D 52, 7182-7204 (1995) doi:10.1103/PhysRevD.52.7182 [arXiv:hep-ph/9506475 [hep-ph]].


# An Improved A* Search Algorithm for Road Networks Using New Heuristic
Estimation

Kevin Y. Chen111Academy for Mathematics, Science, and Engineering, Rockaway,
NJ 07866, USA

###### Abstract

Finding the shortest path between two points in a graph is a fundamental
problem that has been well-studied over the past several decades. Shortest
path algorithms are commonly applied to modern navigation systems, so our
study aims to improve the efficiency of an existing algorithm on large-scale
Euclidean networks. The current literature lacks a deep understanding of
certain algorithms’ performance on these types of networks. Therefore, we
incorporate a new heuristic function, called the $k$-step look-ahead, into the
A* search algorithm and conduct a computational experiment to evaluate and
compare the results on road networks of varying sizes. Our main findings are
that this new heuristic yields a significant improvement in runtime,
particularly for larger networks when compared to standard A*, as well as that
a higher value of $k$ is needed to achieve optimal efficiency as network size
increases. Future research can build upon this work by implementing a program
that automatically chooses an optimal $k$ value given an input network. The
results of this study can be applied to GPS routing technologies or other
navigation devices to speed up the time needed to find the shortest path from
an origin to a destination, an essential objective in daily life.

##  1 Introduction

Finding the shortest path from an origin to a destination is a fundamental,
well-studied problem faced by researchers and mathematicians for at least the
past century. Even in early societies, finding shortest paths was a common,
important objective. However, extensive research in the shortest path problem
began relatively late as compared to other combinatorial optimization
problems, likely due to its apparent simplicity [1]. In recent years, these
algorithms have been applied to various modern navigation systems such as GPS
routing technologies on road networks.

The problem is typically represented by a graph consisting of nodes and edges,
where the goal is to find the shortest path between two given nodes such that
the sum of the weights of the traversed edges is minimized. When applying this
theoretical model to road networks, nodes often represent intersections, and
edges correspond to road segments that are weighted by the length of the
segment. Perhaps the most well-documented shortest path algorithm was
presented by Dutch computer scientist Edsger W. Dijkstra in 1959 [2], which
remains one of the most efficient methods to this day [3]. Dijkstra’s
algorithm finds the optimal solution to the shortest path problem with
nonnegative edge weights by expanding neighboring nodes uniformly in all
directions. Another efficient but less-studied technique is the A* (pronounced
A-star) algorithm [4], which builds off of Dijkstra’s by incorporating a
heuristic function that takes into consideration the estimated remaining
distance from the current node to the destination, thus limiting the search
area.

The heuristic function of the A* search algorithm is a critical part of its
implementation because it highly influences the algorithm’s runtime and
solution optimality [5]. While shortest path problems in general are a well-
studied topic, the current literature notably lacks a detailed analysis on how
the heuristic estimating method of A* affects its computational performance on
real-life networks of various sizes.

Therefore, the primary objective of this work is to introduce a new heuristic
for A* that has not before been implemented and analyzed in the environment of
real road networks. This new look-ahead heuristic involves the idea of
deciding which node in the graph to expand next based on the likelihood that
one of the node’s children is on the shortest path. We will then evaluate the
computational efficiency of this new version of A* and compare it to the
standard A* and Dijkstra’s algorithm by testing them on road networks of three
different sizes. The results of this experiment can be used in the future to
determine which of these algorithms is the most optimal to use on Euclidean
networks under different circumstances.

##  2 Literature Review

While a plethora of past studies have dealt with efforts to compare shortest
path algorithms, many of them evaluate their computational performance on
randomly generated graphs as opposed to real networks [6, 7, 8]. This is
problematic when it comes to real-life applications because synthetic graphs
lack geographical reality and the most efficient shortest path algorithms on
random networks may not practically be ideal when used on empirical systems
[3]. Furthermore, these studies exclude the A* algorithm in their comparisons,
but this decision is somewhat reasonable because A* functions best in the
context of Euclidean space instead of randomized networks. While Zhan and Noon
[9] examine the performance of multiple shortest path algorithms on several
large-scale U.S. road networks, they do not mention A*. Additionally, Chan et
al. [10] compare the efficiency of six different algorithms including A*, but
do so on a relatively small bus route that is not representative of a typical
large transportation network.

On the other hand, other previous experiments carry out direct comparisons
between A* and Dijkstra’s algorithm on complex road networks, but do so
without varying certain features of the algorithms [11, 12]. They generally
find the standard A* to discover the shortest path faster in those
environments. Zeng and Church [3] analyze the effects of varying the type of
data structure used to maintain the set of candidate nodes (e.g. $k$-array
heaps, double buckets, approximate buckets) in A* and Dijkstra’s. Moreover, in
the previously mentioned papers that study the application of shortest path
algorithms to road networks, Euclidean distance is the most commonly used
heuristic function in A*, as it ideally describes the lower bound of the true
cost from a node to the destination in Euclidean space. Some more studies
attempt to incorporate an “overdo” heuristic into A* that sacrifices solution
optimality for a shorter runtime [13, 14]. Finally, others explore a variation
of A* known as the ALT algorithm, which uses a different lower-bounding
technique based on landmarks and the triangle inequality as opposed to
Euclidean bounds [15, 16].

To fill in the aforementioned gaps in the literature, our main contribution is
to improve A*’s performance on large road networks by integrating a heuristic
function that has the potential to decrease runtime while guaranteeing the
optimal path.

##  3 Purpose

The goal of this study is to improve upon the performance of the A* search
algorithm on Euclidean road networks by implementing a $k$-step look-ahead
heuristic. Specifically, we break this problem down into the following steps:

  1. 1.

We aim to evaluate the computational performance of Dijkstra’s algorithm, the
standard A* algorithm, and the A* algorithm with a $k$-step look-ahead
heuristic function on three real road networks of different sizes. We compare
the efficiency of these algorithms in terms of program runtime.

  2. 2.

We will then use the results of the computational experiment to determine
which of the algorithms is best to use for different sizes of the shortest
path problem on Euclidean networks.

##  4 Shortest Path Algorithms

In the past, scientists and mathematicians have derived various algorithms for
solving the shortest path problem. In this paper, we focus on evaluating
Dijkstra’s algorithm and the A* algorithm, and we provide an overview of each
of these methods in the following sections.

###  4.1 Notation and Definitions

In this section, we discuss the notation necessary for understanding graphs
and shortest path algorithms. Most of the notation introduced here will be
used throughout section 4. We then rigorously define the shortest path
problem.

Let $G=(V,E)$ be a directed, undirected, or mixed graph consisting of a set of
vertices (i.e. nodes) $V$ and a set of edges $E$. We denote $e_{x,y}\in E$ as
the edge from node $x$ to node $y$, and $c(x,y)$ as the cost (i.e. weight) of
that edge. We also say that a node $n^{\prime}$ is adjacent to node $n$ if
there exists an edge $e_{n,n^{\prime}}$. In the context of shortest path
algorithms, we let $d(n)$ be the total cost of the shortest path from the
source node to node $n$. We also call this value the distance value of node
$n$. We let the source node of the shortest path problem be $s$ and the target
node be $t$. We also have $\textsf{dist}(x,y)$ represent the Euclidean (i.e.
straight-line) distance between nodes $x$ and $y$.

Next, let

$P=((n_{1},n_{2}),(n_{2},n_{3}),\dots,(n_{m-1},n_{m}))$

be a path of edges connecting node $n_{1}$ to $n_{m}$ where each
$(n_{i},n_{i+1})$ is an edge of the path connecting adjacent nodes $n_{i}$ and
$n_{i+1}$ for all $1\leq i<m$. We say that a path $P$ is valid if and only if
$y_{j}=x_{j+1}$ for all $1\leq j<m-1$, where $(x_{j},y_{j})$ and
$(x_{j+1},y_{j+1})$ are the $j$th and $(j+1)$th elements of $P$, respectively.
We define a path to be acyclic if the path does not pass through any same node
more than once. We call node $n^{\prime}$ a successor of $n$ if $n^{\prime}$
comes after $n$ in a given path. In addition, a path between node $n$ and
$n^{\prime}$ is said to have length $k$ if the path passes through exactly $k$
edges (note that this is a different quantity from total cost of the path).

The goal of the single-pair shortest path problem is to find the path $P$
between $s$ and $t$, where $n_{1}=s$ and $n_{m}=t$, that minimizes the sum
$\sum_{i=1}^{m-1}{c(n_{i},n_{i+1})}$ over all possible values of $m$. Other
common variations of the shortest path problem exist besides the single-pair
problem. These include the single-source problem, where the objective is to
find the shortest paths from one node to all other nodes, and the all-pairs
problem, where a shortest path must be found between all pairs of nodes in the
graph [17].

###  4.2 Dijkstra’s Algorithm

First introduced by Dijkstra in 1959 [2], this algorithm finds the optimal
path on a graph with nonnegative weights from a single source node to a target
or set of targets. While the algorithm is more commonly used to solve the
single-source shortest path problem, it will be implemented in this study to
handle a single pair of nodes for the purpose of comparing its computational
efficiency with that of A*.

Dijkstra’s method, as outlined in Algorithm 1, utilizes a min-priority queue
to store the set of unvisited nodes, leading to more efficient computations
compared to other basic data structures [3]. Since our work focuses on finding
the shortest path from an origin to a particular destination, the algorithm
terminates when the target node $t$ is marked as visited. Asymptotically, the
algorithm has a worst-case time complexity of $O((|V|+|E|)\log|V|)$ when a
standard binary heap is implemented as the priority queue data structure [18].

Algorithm 1 Dijkstra’s Algorithm

$\textsf{Dijkstra}(G,s,t)$:
Input: A graph $G=(V,E)$, start node $s\in V$ and target node $t\in V$.
Output: The shortest path of edges from node $s$ to node $t$.

1:$\textsf{unvisitedQueue}\coloneqq$ priority queue of unvisited nodes sorted
by $d$ value

2:$\textsf{previous}\coloneqq$ an empty map

3:$d(n)\coloneqq\infty$ for all nodes $n$ in $G$

4:$d(s)\coloneqq 0$

5:while unvisitedQueue not empty and $t$ not visited do

6: Extract node $i$ from unvisitedQueue

7: for all $j$ adjacent to $i$ do

8: if $j$ is visited then

9: Continue

10: $\textsf{dTemp}\coloneqq d(i)+c(i,j)$

11: if $\textsf{dTemp}<d(j)$ then

12: $d(j)\coloneqq\textsf{dTemp}$

13: $\textsf{previous}[j]\coloneqq i$

14: Reorder unvisitedQueue

###  4.3 A* Algorithm

The A* search algorithm is a best-first search method [19] and has a general
procedure similar to Dijkstra’s, but the difference comes with the concept of
a heuristic function that estimates the remaining distance from node $n$ to
the target node $t$. This crucial feature causes A* to run faster than
Dijkstra’s in many past computational experiments [3, 11, 13]. A* functions by
choosing which path to extend based on the node that minimizes the value of

$f(n)=g(n)+h(n),$

where $n$ is the next node in the path, $g(n)$ is the current distance value
of $n$, and $h(n)$ is the predetermined heuristic value of $n$. Using this
heuristic, A* is able to concentrate its search in the direction of the target
by prioritizing nodes with a lower $h$ value (i.e. nodes that are estimated to
be closer to the destination). This main distinction between A* and Dijkstra’s
can be seen in Algorithm 2.

Algorithm 2 A* Search Algorithm

$\textsf{A*}(G,s,t,h)$:
Input: A graph $G=(V,E)$, start node $s\in V$, target node $t\in V$, and
heuristic function $h(n)$.
Output: The shortest path of edges from node $s$ to node $t$.

1:$\textsf{unvisitedQueue}\coloneqq$ priority queue of unvisited nodes sorted
by $f$ value

2:$\textsf{previous}\coloneqq$ an empty map

3:$f(n)\coloneqq\infty$ for all nodes $n$ in $G$

4:$g(n)\coloneqq\infty$ for all nodes $n$ in $G$

5:$f(s)\coloneqq h(s)$

6:$g(s)\coloneqq 0$

7:while unvisitedQueue not empty and $t$ not visited do

8: Extract node $i$ from unvisitedQueue

9: for all $j$ adjacent to $i$ do

10: if $j$ is visited then

11: Continue

12: $\textsf{gTemp}\coloneqq g(i)+c(i,j)$

13: if $\textsf{gTemp}<g(j)$ then

14: $g(j)\coloneqq\textsf{gTemp}$

15: $f(j)\coloneqq g(j)+h(j)$

16: $\textsf{previous}[j]\coloneqq i$

17: Reorder unvisitedQueue

The most commonly used heuristic estimator when dealing with real networks is
Euclidean distance, as this metric always describes the lowest bound of the
distance between any two points in space [20, 21]. Euclidean distance is an
example of an admissible heuristic, which is defined as a heuristic that never
overestimates the true lowest possible cost from a node $n$ to $t$, thereby
guaranteeing that the shortest path to $t$ is not overlooked. Additionally, a
heuristic is considered consistent (or monotone) if it satisfies

$h(n)\leq c(n,p)+h(p)\text{ and }h(t)=0,$

where $p$ is a successor of any node $n$ in the graph. In the context of A*,
this means that the cost by which a node is reached must be the lowest
possible once that node is expanded, given that the graph only includes
nonnegative edge weights. When the heuristic is simply Euclidean distance,
this property of consistency is analogous to the triangle inequality because
it essentially states that the straight-line distance between $n$ and $t$ is
never greater than the sum of the distance between $n$ and $p$ plus the
distance between $p$ and $t$. A consistent heuristic is always admissible, but
the converse is not necessarily true. This can be proven by induction [22].

While the time complexity of A* depends heavily on the search space and the
heuristic function used, its worst-case performance on an unbounded search
space is $O(b^{d})$, where $b$ is the branching factor (i.e. the average
number of children per node, or the outdegree) and $d$ is the depth of the
shortest path solution [19].

###  4.4 The $k$-step Look-ahead Heuristic

We aim to improve the computational efficiency of the A* search algorithm by
integrating the $k$-step look-ahead heuristic, in which $k$ is a parameter
that we choose. First, we start with a couple of definitions. We define
$C(P,t)$ to be a cut function that represents the section of path $P$ up until
node $t$, if $t$ is on path $P$. For instance, if
$P=((a,b),(b,t),(t,c),(c,d))$, then $C(P,t)=((a,b),(b,t))$. If $P$ does not
pass through $t$, then $C(P,t)=P$. If $X$ is the set of all visited nodes in
$G$, then let $S_{G,X,n,k}$ be the set of all acyclic paths from $n$ of length
$k$ that do not pass through any nodes in $X$.

The intuitive explanation behind the $k$-step look-ahead heuristic is that the
heuristic scans next possible paths to a certain extent from a given node and
selects the estimated cost of the path that seems most promising.
Specifically, instead of estimating the remaining cost to get from node $n$ to
target node $t$ directly using Euclidean distance, we do so by iteratively
finding the cost of each possible path of length $k$ to the unvisited
successors of $n$, added to the Euclidean distance between that successor node
and $t$. We then select the lowest of those calculated values and assign that
to $h(n)$. More formally, we have

$h(n)=\min_{P\in S_{G,X,n,k}}\left[\left(\sum\limits_{(x,y)\in
C(P,t)}c(x,y)\right)+\textsf{dist}(n^{\prime},t)\right]$

where $n^{\prime}$ is the final node in path $C(P,t)$.

For example, to find the estimated cost from $n$ to $t$ using the 1-step look-
ahead heuristic, we calculate

$c(n,n^{\prime})+\textsf{dist}(n^{\prime},t)$

for every unvisited child node $n^{\prime}$ of $n$, and assign $h(n)$ to the
smallest resulting value.

We emphasize two details implied from our definition of the $k$-step look-
ahead heuristic that become increasingly important as the value of $k$ gets
larger:

  1. 1.

As previously mentioned, the “look-ahead” path $P$ that our algorithm takes
must be acyclic, i.e. it should not revisit any nodes previously passed
through earlier in the same look-ahead path. Otherwise, the process could be
simulating a sub-optimal path where nodes are revisited, potentially leading
to unnecessary calculations. This is particularly true when dealing with
undirected graphs where it is possible to travel between nodes in both
directions. The list alreadySeen in Algorithm 3 keeps track of which nodes
have already been looked through in a given look-ahead path.

  2. 2.

If the algorithm reaches the target node during its look-ahead process, that
particular look-ahead path should terminate and the returned value should be
the cost of that path from node $n$ to $t$. This prevents the algorithm from
unnecessarily extending to successor nodes beyond the desired target node.
This is the reasoning behind finding the cost of path $C(P,t)$ as opposed to
the entire cost of $P$.

Algorithm 3 shows example pseudocode for computing the 2-step look-ahead
heuristic value of a node $n$.

Algorithm 3 2-step Look-Ahead Heuristic

$\textsf{2-step Look-ahead}(G,n,t)$:
Input: A Euclidean graph $G=(V,E)$, node $n\in V$, and target node $t\in V$.
Output: $h(n)$, or the estimated cost of the shortest path from $n$ to $t$.

1:if $n=t$ then

2: $h(n)\coloneqq 0$

3:else

4: $\textsf{alreadySeen}\coloneqq$ an empty list only containing $n$

5: $\textsf{possibleVals}\coloneqq$ an empty list only containing $\infty$

6: for all $n_{1}$ adjacent to $n$ do

7: if $n_{1}$ has not been visited then

8: if $n_{1}=t$ then

9: Add $c(n,n_{1})$ to possibleVals

10: Continue

11: Add $n_{1}$ to alreadySeen

12: for all $n_{2}$ adjacent to $n_{1}$ do

13: if $n_{2}$ is not visited and not in alreadySeen then

14: Add $c(n,n_{1})+c(n_{1},n_{2})+\textsf{dist}(n_{2},t)$ to possibleVals

15: Remove $n_{1}$ from alreadySeen

16: $h(n)\coloneqq\min$(possibleVals)

17:return $h(n)$

The $k$-step look-ahead is a more accurate heuristic of the shortest path
between $n$ and $t$ than direct Euclidean distance because it takes into
account some of the edges that the real shortest path from $s$ to $t$ will
actually traverse. This suggests that a shortest path will be found quicker
with this new heuristic than by using standard A*, particularly for relatively
small values of $k$ (for larger $k$, the number of calculations the heuristic
has to make tends to increase exponentially, which can lead to slower
runtimes). Furthermore, this estimator is still admissible, as the resulting
value for $h(n)$ will never be greater than the cost of the true shortest path
from $n$ to $t$, so the optimal solution is guaranteed on Euclidean networks.

##  5 Methodology

In this section, we discuss the procedure used to perform our computational
experiment on the efficiency of Dijkstra’s algorithm, standard A*, and A* with
a $k$-step look-ahead heuristic for various values of $k$ on sample road
networks. To our knowledge, the new heuristic has not before been evaluated on
Euclidean networks. All tests are run using the Python programming language
[23] in Google Colaboratory on a computer with a 1.3 GHz Intel Core i7
processor and 16 GB of RAM.

###  5.1 Test Network Details

We implement the shortest path algorithms on three different road networks as
described in Table 1: the relatively small Downtown Brooklyn network, the
medium-sized Jersey City network, and the larger New York City network (see
Appendix A for visual layouts). We choose to include datasets of different
sizes for the purpose of determining which algorithms perform comparatively
well or poorly for various scales of networks. The data is obtained from the
OSMnx Street Networks Dataverse, a compilation found on Harvard Dataverse that
includes the street networks of every city, town, urban area, county, census
tract, and Zillow-defined neighborhood in the United States [24, 25]. The
networks are created from OpenStreetMap data using the OSMnx software [26]. In
the networks, nodes correspond to street intersections, edges represent the
road segments connecting those intersections, and weights correspond to the
lengths of those road segments in meters. Additionally, for each network, the
data source provides node and edge lists. The node lists contain each node in
the network, its assigned ID, and its geographic coordinates, while the edge
lists include information about the two nodes that each edge connects, the
cost of the edge, and whether the edge is one-way only, among other features.
Note that these networks are considered mixed graphs because they contain both
directed and undirected edges.

Table 1: Characteristics of test road networks. Road network | # of nodes | # of edges
---|---|---
Downtown Brooklyn | 149 | 281
Jersey City | 2232 | 5320
New York City | 55136 | 141085

Each road network dataset must also be pre-processed in order to remove groups
of nodes that are entirely disconnected from the rest of the network, thus
ensuring a fully navigable road system. Lastly, before the final experiment is
run, all algorithms are tested on several source-target pairs for all three
networks to confirm that they consistently find the same optimal path.

###  5.2 Algorithm Implementations

In our computational experiment, twelve algorithms are tested on each network:
Dijkstra’s, standard A*, and A* with a $k$-step look-ahead heuristic where $k$
is each of 1 to 10 inclusive. For each network, a certain number of source-
target node pairs are randomly selected from the network (with the exception
of New York City), and each algorithm is run a fixed number of times for each
source-target pair. Specifics are outlined in Table 2. We then record the
resulting average runtime and average nodes expanded (i.e. nodes visited) for
an algorithm’s execution (see Results section). As Table 2 shows, we design
our experiment differently between road networks in terms of the number of
times each algorithm is run. This is because the overall runtimes of the
algorithms significantly increase with the size of the network, and runtimes
on the smaller networks must be recorded at a greater degree of precision.
Since we run the algorithms on only one distinct shortest path problem in the
New York City dataset, we deliberately choose the source and destination nodes
to be on opposite sides of the network in order to encapsulate the network’s
size and complexity.

Table 2: Experimental details for each test road network.

Road Network | # of source-target pairs | # of runs per pair
---|---|---
Downtown Brooklyn | 200 | 30
Jersey City | 100 | 10
New York City | 1 | 1

To implement a priority queue data structure to efficiently store unvisited
nodes, we utilize Python’s heapq module, which represents a binary tree and
can order a list by priority in linear time. Moreover, we incorporate a
calculation of Euclidean distance between two given nodes in meters using the
haversine formula since the geographic coordinates of each node is known [12].
The formula states that

$d=2r\arcsin\left(\sqrt{\sin^{2}\left(\dfrac{\phi_{2}-\phi_{1}}{2}\right)+\cos(\phi_{1})\cos(\phi_{2})\sin^{2}\left(\dfrac{\lambda_{2}-\lambda_{1}}{2}\right)}\right),$

where

  * •

$(\phi_{1},\lambda_{1})$ and $(\phi_{2},\lambda_{2})$ are the coordinates of
the two points in radians,

  * •

$r$ is the radius of the Earth, and

  * •

$d$ is the great-circle distance between the two points.

Finally, for our implementations of all variations of A*, the final recorded
runtimes do not include the program’s graph pre-processing steps of
calculating the Euclidean distance between each node in the network and the
target node.

##  6 Results

The results of our computational experiment with shortest path algorithms on
road networks are displayed in this section. As previously described, we run
the A* algorithm with a $k$-step look-ahead heuristic on each of the three
sample road networks for all $k$ from 1 to 10 inclusive. Additionally, we
include Dijkstra’s algorithm to provide context of the algorithms’
performance. As is a common convention [3, 9, 10, 11], we use runtime to judge
an algorithm’s efficiency. For each source-target pair, all algorithms we
implement guarantee the same shortest path. See Appendix B for additional
visualizations of the number of nodes expanded as a function of $k$.

###  6.1 Small Sample Network

Results for all twelve algorithms on the Downtown Brooklyn street network are
shown in Table 3. In terms of runtime, a 2-step look-ahead heuristic achieves
the greatest efficiency by a small margin for this particular network at
roughly 0.0026 seconds on average and 10.5% faster than standard A*. Moreover,
the heuristics with $k=1$ and $k=3$ also perform better than regular A*.
Figure 1 shows that as $k$ increases, A*’s average runtime grows
exponentially. While $k=2$ achieves the quickest runtime, A* expands the least
number of nodes on average with a 4-step look-ahead. In addition, when $k=10$,
the largest value of $k$ that we test, the runtime is significantly greater
than that of Dijkstra’s while still expanding less than half as many nodes.

Table 3: Comparison of Algorithm Performance on Downtown Brooklyn Network.

Algorithm | Avg Runtime (s) | Avg # of Nodes Expanded
---|---|---
Dijkstra’s | 0.0055 | 69.81
Standard A* | 0.0030 | 34.44
A*, $k=1$ | 0.0028 | 30.59
A*, $k=2$ | 0.0026 | 27.83
A*, $k=3$ | 0.0029 | 26.92
A*, $k=4$ | 0.0032 | 25.59
A*, $k=5$ | 0.0048 | 26.91
A*, $k=6$ | 0.0057 | 28.98
A*, $k=7$ | 0.0089 | 29.39
A*, $k=8$ | 0.0119 | 29.34
A*, $k=9$ | 0.0225 | 31.28
A*, $k=10$ | 0.0368 | 30.68

Figure 1: A visualization of the average algorithm runtime in seconds as a
function of $k$ in A* implementation with a $k$-step look-ahead heuristic when
run on the Downtown Brooklyn network. A value of $k=0$ indicates the standard
A* algorithm. Optimal efficiency for this network is achieved with a 2-step
look-ahead, but its improvement from standard A* is negligible.

###  6.2 Medium Sample Network

As shown in Table 4, the A* algorithm efficiency is able to be improved for 6
different values of $k$ on the medium-sized Jersey City road network. Most
notably, the best performance occurs when $k=4$, running 32.6% faster than
regular A* at approximately 0.2719 seconds on average. This same value of $k$
also results in the least average number of nodes expanded. We also see a
growing exponential behavior in algorithm runtime for values of $k$ greater
than $4$, depicted in Figure 2. Once $k$ hits 10, A*’s runtime exceeds that of
Dijkstra’s even though it expands much fewer nodes.

Table 4: Comparison of Algorithm Performance on Jersey City Network.

Algorithm | Avg Runtime (s) | Avg # of Nodes Expanded
---|---|---
Dijkstra’s | 1.1554 | 1192.50
Standard A* | 0.4033 | 304.88
A*, $k=1$ | 0.3172 | 247.29
A*, $k=2$ | 0.3196 | 220.90
A*, $k=3$ | 0.2759 | 203.58
A*, $k=4$ | 0.2719 | 192.81
A*, $k=5$ | 0.3088 | 210.99
A*, $k=6$ | 0.3547 | 220.63
A*, $k=7$ | 0.4190 | 233.78
A*, $k=8$ | 0.7356 | 240.70
A*, $k=9$ | 0.9383 | 254.02
A*, $k=10$ | 1.5565 | 261.78

Figure 2: Algorithm runtime in seconds as a function of $k$ in our
implementation of A* search algorithm with a $k$-step look-ahead heuristic
when run on the Jersey City network. A* achieves optimal efficiency using a
4-step look-ahead heuristic for this particular network.

###  6.3 Large Sample Network

The street system of New York City is by far the largest of the three sample
networks we use in this experiment, and all runtimes surpass 10 minutes. We
test each algorithm once on a chosen source-target pair in this network, and
find that an 8-step look-ahead heuristic achieves optimal performance at just
684.10 seconds, as indicated in Table 5. Using this heuristic results in a
22.2% shorter runtime than using standard A*. In total, 9 different values of
$k$ perform better than regular A*. Figure 3 visualizes this data, showing
that the runtimes stay relatively low when $k$ is between 3 and 8, but
skyrockets once $k$ exceeds 8. Additionally, A* expands the least number of
nodes in the New York City network when using a 4-step look-ahead heuristic.
When $k=10$, A* runs more efficiently than Dijkstra’s while also expanding
less than half as many nodes.

Table 5: Comparison of Algorithm Performance on New York City Network.

Algorithm | Runtime (s) | # of Nodes Expanded
---|---|---
Dijkstra’s | 1389.98 | 54346
Standard A* | 879.43 | 23370
A*, $k=1$ | 829.41 | 20498
A*, $k=2$ | 771.66 | 19200
A*, $k=3$ | 722.54 | 18323
A*, $k=4$ | 760.31 | 17769
A*, $k=5$ | 712.97 | 17934
A*, $k=6$ | 811.00 | 18959
A*, $k=7$ | 761.65 | 20418
A*, $k=8$ | 684.10 | 20559
A*, $k=9$ | 850.23 | 20264
A*, $k=10$ | 1047.57 | 20348

Figure 3: Algorithm runtime in seconds as a function of $k$ for A* with a
$k$-step look-ahead heuristic when run on the New York City network. A*
achieves optimal efficiency using an 8-step look-ahead heuristic for this
particular network, but also performs relatively well when $k=3$ and $k=5$.

##  7 Discussion

Through computational experimentation, our goal has been to determine whether
implementing a $k$-step look-ahead heuristic would improve the efficiency of
the A* search algorithm and to decide which algorithms work best on different
sizes of road network. Here, we address this problem by interpreting the
results from the previous section.

First, the results agree with our hypothesis that the $k$-step look-ahead
heuristic would lead to shorter runtime than using direct Euclidean distance
as the estimator. For the Downtown Brooklyn network, although we are able to
slightly improve efficiency for a few values of $k$, the differences in
runtime are negligible ten-thousandths of a second, which is unlikely to be
impactful in the real world. On the other hand, using the new heuristic on the
Jersey City and New York City networks shortens runtime by at most 0.13 and
195.33 seconds, respectively. In addition, on the Downtown Brooklyn network, 3
values of $k$ improve the efficiency of A*, whereas this number is 6 for the
Jersey City network and 9 for New York City. The optimal efficiency occurs
when $k=2$ for Downtown Brooklyn, $k=4$ for Jersey City, and $k=8$ for New
York City. This suggests not only that the $k$-step look-ahead heuristic is
more effective for larger networks, but that the single optimal value of $k$
increases as the network size grows. The above conclusions seem reasonable
because in a larger graph, there is a greater need to expand an algorithm’s
search in the direction of the target node, which is accomplished by greater
values of $k$ to create a more accurate heuristic.

A reason that one shortest path algorithm may run faster than another may be
that it visits and expands less nodes, but this metric does not tell the whole
story in all cases. For example, in each of the two smaller networks (i.e.
Downtown Brooklyn and Jersey City), A* with a 10-step look-ahead heuristic has
a considerably longer runtime than Dijkstra’s algorithm, but the number of
nodes that A* expands in these cases is significantly lower than that of
Dijkstra’s. This is explained by the sheer number of steps that A* has to take
when calculating the cost of each look-ahead path. As mentioned earlier, the
number of steps needed to calculate a $k$-step look-ahead heuristic value
grows exponentially as $k$ increases due to the many possible look-ahead paths
that diverge at each node.

Meanwhile, the situation is different with the New York City network. In this
case, both the runtime and number of nodes expanded are less for A* with a
10-step look-ahead than for Dijkstra’s. This further emphasizes A*’s
suitability for large-scale networks, even for relatively high values of $k$.

A more surprising result of this experiment emerges when looking at data for
the number of nodes expanded. The original speculation is that a greater value
of $k$ leads to a more accurate heuristic, thus expanding fewer “incorrect”
nodes that are not on the shortest path. However, this does not seem to be the
case. In particular, for all three networks, the minimal number of nodes
expanded during A* search occurs when $k=4$. Finding a justification to
explain this phenomenon can be one direction of future work. Another common
theme between all three test networks is that all variations of A* expand
significantly fewer nodes than Dijkstra’s does, likely due to Dijkstra’s
tendency to blindly expand nodes in all directions until the destination is
found, as it is not guided by a heuristic. For instance, in the New York City
network, Dijkstra’s algorithm expands 54346 out of 55136 nodes in the entire
graph, while all A* implementations go about finding the shortest path more
intelligently. Note that the reason Dijkstra’s expands nearly all the nodes in
the network for this specific case is that we deliberately choose the source
and target nodes to be on opposite ends of the network.

While many useful insights are gained from our interpretations, these are
preliminary results that may need more testing to solidify. For example, our
single run on the New York City dataset for each algorithm may not be fully
representative of all possibilities on the network; nevertheless, our analysis
in this section is likely safe to generalize and apply to the real world. To
further strengthen the validity of our conclusions, future experiments should
run these algorithms on many more source-target pairs on several different
road networks.

Overall, the results of this study show that incorporating a $k$-step look-
ahead heuristic into the A* search algorithm yields greater efficiency in
solving the shortest path problem on transportation networks and can be
particularly effective on larger, more complex road systems.

##  8 Conclusion and Future Work

The shortest path problem is a fundamental combinatorial optimization problem
that has been intensely studied for the past several decades. In the current
literature, few studies have compared various shortest path algorithms’
performance on large, complex road networks [3, 9, 13]. The main contribution
of our work is to implement the A* search algorithm with a $k$-step look-ahead
heuristic to solve the single-pair shortest path problem on various road
networks for the first time. We implement Dijkstra’s algorithm, A*, and A*
with a $k$-step look-ahead heuristic for all $k$ from 1 to 10 inclusive. We
then compare these algorithms’ performance on three real road networks of
different sizes: the small Downtown Brooklyn network, the medium-sized Jersey
City network, and the larger New York City network. From the results of our
experiment, we conclude that using the look-ahead heuristic significantly
improves the efficiency of A*, especially on larger networks. We also find
that the greater the network size, the higher $k$ needs to be in order to
achieve optimal efficiency on the given network. Limitations of our
experimental methodology include the lack of a diverse range of source-target
pairs to fully represent all possible shortest path problems on the New York
City network. The results of this study can be incorporated into cutting-edge
navigation systems on transportation networks, such as GPS routing
technologies, as finding the shortest possible path between two locations is a
common everyday objective.

Some possibilities of related future research directions include incorporating
different data structures (i.e. double buckets, approximate buckets) for the
priority queue in the implementation of A* with $k$-step look-ahead, as
opposed to the simple binary heap structure we use in this study.
Additionally, the average degree of a graph (i.e. the average number of edges
that are incident to each node in the graph) may be another factor that
affects the efficiency of using a $k$-step look-ahead heuristic, and should be
analyzed in the context of real road networks. Finally, future studies may
explore the possibility of implementing a program that automatically chooses
the most optimal value of $k$ to use in A*. The efficiency of this algorithm
can then be evaluated on various Euclidean networks.

##  9 Acknowledgements

The author thanks Jesse Stern (Department of Computer Science, University of
Chicago) for feedback and mentorship on this project. The author also thanks
the Summer STEM Institute (SSI).

## References

  * [1] A. Schrijver, “On the history of the shortest path problem,” Documenta Mathematica, vol. 17, no. 1, pp. 155–167, 2012.
  * [2] E. W. Dijkstra and others, “A note on two problems in connexion with graphs,” Numerische mathematik, vol. 1, no. 1, pp. 269–271, 1959.
  * [3] W. Zeng and R. L. Church, “Finding shortest paths on real road networks: the case for A,” International journal of geographical information science, vol. 23, no. 4, pp. 531–543, 2009.  Publisher: Taylor & Francis.
  * [4] P. E. Hart, N. J. Nilsson, and B. Raphael, “A formal basis for the heuristic determination of minimum cost paths,” IEEE transactions on Systems Science and Cybernetics, vol. 4, no. 2, pp. 100–107, 1968.  Publisher: IEEE.
  * [5] L. Fu, D. Sun, and L. R. Rilett, “Heuristic shortest path algorithms for transportation applications: State of the art,” Computers & Operations Research, vol. 33, no. 11, pp. 3324–3343, 2006.  Publisher: Elsevier.
  * [6] B. V. Cherkassky, A. V. Goldberg, and T. Radzik, “Shortest paths algorithms: Theory and experimental evaluation,” Mathematical programming, vol. 73, no. 2, pp. 129–174, 1996.  Publisher: Springer.
  * [7] P. Biswas, P. Mishra, and N. Mahanti, “Computational Efficiency of Optimized Shortest Path Algorithms.,” Int. J. Comput. Sci. Appl., vol. 2, no. 2, pp. 22–37, 2005.  Publisher: Citeseer.
  * [8] K. Magzhan and H. M. Jani, “A review and evaluations of shortest path algorithms,” International journal of scientific & technology research, vol. 2, no. 6, pp. 99–104, 2013.
  * [9] F. B. Zhan and C. E. Noon, “Shortest path algorithms: an evaluation using real road networks,” Transportation science, vol. 32, no. 1, pp. 65–73, 1998\.  Publisher: INFORMS.
  * [10] S. Y. M. Chan, N. Adnan, S. S. Sukri, and W. M. N. Wan Zainon, “An experiment on the performance of shortest path algorithm,” in Knowledge Management International Conference (KMICe), pp. 7–12, 2016.
  * [11] D. Rachmawati and L. Gustin, “Analysis of Dijkstra’s Algorithm and A* Algorithm in Shortest Path Problem,” in Journal of Physics: Conference Series, vol. 1566, p. 012061, IOP Publishing, 2020.  Issue: 1.
  * [12] J. Jain, “Shortest Path for Emergency Services: A Comparative Analysis,” International Journal for Research in Applied Science and Engineering Technology, vol. 9, pp. 426–430, June 2021.
  * [13] F. Engineer, “Fast shortest path algorithms for large road networks,” in Proceedings of 36th annual ORSNZ conference, 2001.
  * [14] R. Jacob, M. Marathe, and K. Nagel, “A computational study of routing algorithms for realistic transportation networks,” Journal of Experimental Algorithmics (JEA), vol. 4, pp. 6–es, 1999.  Publisher: ACM New York, NY, USA.
  * [15] A. V. Goldberg and C. Harrelson, “Computing the shortest path: A search meets graph theory.,” in SODA, vol. 5, pp. 156–165, Citeseer, 2005.
  * [16] A. V. Goldberg, H. Kaplan, and R. F. Werneck, “Reach for A*: Shortest Path Algorithms with Preprocessing.,” in The shortest path problem, pp. 93–139, Citeseer, 2006.
  * [17] A. Madkour, W. G. Aref, F. U. Rehman, M. A. Rahman, and S. Basalamah, “A survey of shortest-path algorithms,” arXiv preprint arXiv:1705.02044, 2017\.
  * [18] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, Introduction to algorithms.  MIT press, 3rd ed., 2009.
  * [19] S. Russell and P. Norvig, Artificial Intelligence: A Modern Approach.  Pearson Education, 3rd ed., 2010.
  * [20] S. K. Sharma and B. Pal, “Shortest path searching for road network using a* algorithm,” International Journal of Computer Science and Mobile Computing, vol. 4, no. 7, pp. 513–522, 2015.
  * [21] D. Schultes, “Route Planning in Road Networks.,” in Ausgezeichnete Informatikdissertationen, pp. 22–29, 2008.
  * [22] R. Dechter and J. Pearl, “Generalized best-first search strategies and the optimality of A,” Journal of the ACM (JACM), vol. 32, no. 3, pp. 505–536, 1985.  Publisher: ACM New York, NY, USA.
  * [23] K. Hong, “Dijkstra’s shortest path algorithm.” https://www.bogotobogo.com/python/python_Dijkstras_Shortest_Path_Algorithm.php, 2020\.
  * [24] G. Boeing, “U.S. Street Network Shapefiles, Node/Edge Lists, and GraphML Files,” 2017.
  * [25] G. Boeing, “A multi-scale analysis of 27,000 urban street networks: Every us city, town, urbanized area, and zillow neighborhood,” Environment and Planning B: Urban Analytics and City Science, vol. 47, no. 4, pp. 590–608, 2020\.
  * [26] G. Boeing, “Osmnx: New methods for acquiring, constructing, analyzing, and visualizing complex street networks,” Computers, Environment and Urban Systems, vol. 65, pp. 126–139, 2017.

##  A Sample Road Networks

Figure 4: Downtown Brooklyn network (149 nodes, 281 edges). Figure 5: Jersey
City network (2232 nodes, 5320 edges). Figure 6: New York City network (55136
nodes, 141085 edges).

##  B Nodes Expanded vs. $k$ Graphs

Figure 7: Average number of nodes expanded as a function of $k$ on the
Downtown Brooklyn network. Minimal nodes are expanded at $k=4$. Figure 8:
Average number of nodes expanded as a function of $k$ on the Jersey City
network. Minimal nodes are expanded at $k=4$. Figure 9: Average number of
nodes expanded as a function of $k$ on the New York City network. Minimal
nodes are expanded at $k=4$.


* Veilleux et al. (2006) Veilleux, S., Kim, D. C., Peng, C. Y., et al. 2006, ApJ, 643, 707, doi: 10.1086/503188
  * Veilleux et al. (2009) Veilleux, S., Rupke, D. S. N., Kim, D. C., et al. 2009, ApJS, 182, 628, doi: 10.1088/0067-0049/182/2/628
  * Veilleux et al. (2013) Veilleux, S., Meléndez, M., Sturm, E., et al. 2013, ApJ, 776, 27, doi: 10.1088/0004-637X/776/1/27
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wada et al. (2018) Wada, K., Fukushige, R., Izumi, T., & Tomisaka, K. 2018, ApJ, 852, 88, doi: 10.3847/1538-4357/aa9e53
  * Wada et al. (2016) Wada, K., Schartmann, M., & Meijerink, R. 2016, ApJ, 828, L19, doi: 10.3847/2041-8205/828/2/L19
  * Wheeler et al. (2020) Wheeler, J., Glenn, J., Rangwala, N., & Fyhrie, A. 2020, ApJ, 896, 43, doi: 10.3847/1538-4357/ab8f32
  * Winter et al. (2012) Winter, L. M., Veilleux, S., McKernan, B., & Kallman, T. R. 2012, ApJ, 745, 107, doi: 10.1088/0004-637X/745/2/107
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Yamada et al. (2013) Yamada, R., Oyabu, S., Kaneda, H., et al. 2013, PASJ, 65, 103, doi: 10.1093/pasj/65.5.103
  * Yamashita et al. (2017) Yamashita, T., Komugi, S., Matsuhara, H., et al. 2017, ApJ, 844, 96, doi: 10.3847/1538-4357/aa7af1
  * Yang et al. (2010) Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062, doi: 10.1088/0004-637X/718/2/1062
  * Yuan et al. (2010) Yuan, T. T., Kewley, L. J., & Sanders, D. B. 2010, ApJ, 709, 884, doi: 10.1088/0004-637X/709/2/884