# Eccentric Pairs: Analytic Gravitational Waves from Binary Black Holes in
Elliptic Orbits

Dillon Buskirk 111Present address: Department of Physics & Astronomy,
University of Kentucky, Lexington, KY 40506, USA and Maria C. Babiuc Hamilton
222Corresponding author<EMAIL_ADDRESS>Department of Physics, Marshall
University, Huntington, WV 25755, USA

###### Abstract

Gravitational waves (GW) from eccentric binaries have intricate signals
encoding important features about the location, creation and evolution of the
sources. Eccentricity shortens the merger time, making the emitted GW
statistically predominant in the observed data once detectors will reach the
required sensitivity. We present a novel implementation of fully analytical GW
templates from eccentric binary black hole (BBH) mergers within the Wolfram
Mathematica software. We increase the accuracy by identifying and minimizing
the possible source of errors. We start with an overview of the physics
involved in eccentric mergers, then assemble the strain for the inspiral by
employing up to six post-Newtonian (PN) corrections. We complete the eccentric
inspiral with the quasi-circular Backwards one Body (BoB) merger model in
frequency, amplitude and phase, then we build the hybrid GW strain for the
whole evolution of the binary. For low eccentricity we reach coincidence in
the overlap, with no ambiguity in the time interval, a remarkable improvement
from the usual matching techniques. For high-eccentricity we compensate for
the implicit quasi-circular assumption of the BOB approach, by introducing a
small rescaling in amplitude. Our streamlined implementation is relevant for
the new field of GW astronomy and is straightforward to understand, use and
extend, offering researchers in the field a valuable open resource tool.

  *

> If you thought that science was certain - well, that is just an error on
> your part.
>
> Richard P. Feynman

“Not only is the Universe stranger than we think, it is stranger than we can
think.”

##  1 Introduction

In the 17th century, Johannes Kepler introduced elliptical orbits and his
three famous laws of planetary motion. Since then, it is universally accepted
that all binary systems move in elliptical orbits around the center of mass.
This empirical observation was used by Newton to build the law of universal
gravitation, unchallenged for three centuries, until Einstein introduced the
theory of general relativity (GR) in 1915. This theory explained the observed
slight precession in the elliptical orbit of Mercury, known as the
relativistic perihelion advance. More important, it showed that a binary is
bound to spiral inward and circularize by releasing orbital energy as
gravitational waves (GW) until the inevitable collision occurs. It may seem
like a rare scenario that two compact objects lock each other into a
collapsing orbit. However, stars are usually born in pair, and more than
$70\%$ of their population is massive enough to end up either as a black hole
or a neutron star [1].

By the time they reach Earth, GWs have extremely small amplitudes, with an
expected strain on a given length of about $10^{-20}$. This makes their
detection very challenging and possible only for very strong gravity, as in
the case of the two merging black holes first detected by the advanced Laser
Interferometer Gravitational-Wave Observatory (aLIGO) in September 2015 [2].
This discovery created a new branch of science, the Gravitational-wave
Astronomy [3]. The data and software employed in those detections are freely
available [4], enabling peer scientists to test and replicate them [5], an
essential requirement for the upholding of high scientific standards.

All the detections reported until now come from circularized binaries, but
this is not necessarily because they had zero eccentricity [6]. The emission
of GWs decreases the orbital eccentricity [7], thus very small eccentricity is
expected close to the merger, when the signal enters the sensitivity band of
the current detectors, especially for isolated compact binaries evolving for a
long time. However, in certain populations of merging binaries living in dense
stellar environments, such as galactic cores and globular clusters [8, 9], or
in isolated triplets [10], eccentricity is present even at small separations.
Current ground-based detectors such as aLIGO and aVirgo, did not reach yet the
sensitivity levels required for the detection of eccentricity. They are only
tuned to observe the brightest GW sources, in a limited frequency band, very
close to collision [11]. In the near future though, as the technical
difficulties are solved, detection of eccentricities will become possible
[12]. Moreover, next generation of gravitational-wave observatories are being
built and will come online within the next decade [13, 14, 15, 16, 17, 18].
Those instruments will be able to detect GWs with up to 10 times increase in
sensitivity, to discern signals at an earlier stage in the life of the binary,
when orbital eccentricity is still large and to track their evolution for a
long time. It is thus only a matter of time until eccentric GWs will be
detected.

GWs with eccentricity are much richer in structure and they reveal valuable
information on how the binary came together, the medium in which the pair
resides and its long-time evolution [19, 20, 21]. It is known that
eccentricity decreases the merger time, meaning that GWs with eccentricity
encoded in the signal will become ubiquitous once detectors will have the
demanded sensitivity. But to decipher any such detection and to determine
their proprieties, it is essential to have accurate templates that are valid
for small to moderate eccentricities.

GW templates are generated mainly using two families of models currently
implemented in the LIGO Scientific Collaboration Algorithm Library (LALSuite)
[LAL], the Effective-one-Body approach (EOB) [22, 23, 24] and the
phenomenological frequency domain model (IMRPhenom) with its flavors [25, 26,
27, 28]. Those waveform models are following in essence the same two-steps
procedure to construct complete templates for the signal, from the time it
enters the detector’s frequency band, through the peak of the amplitude (or
chirp) emitted when the two black holes collide, until it gradually becomes
smaller, and ending when the final black hole settles down. First, the binary
motion is split in three distinct regions: inspiral, merger, and ringdown,
then different mathematical formalisms are applied to obtain the waveform for
each region. Finally, a complete GW template is generated for the whole
evolution of the binary, by matching those regions in time or in frequency.

The inspiral - or the weak field region - is modeled using different versions
of the post-Newtonian (PN) formalism, which obtains solutions to Einstein’s
equations of GR as higher-order expansions of Newton’s law of gravity in terms
of the small PN parameter $x_{PN}=v^{2}/c^{2}$, also called the _slowness_.
Here, $v$ is the orbital velocity and $c$ the speed of light (see [29, 30, 31,
32, 33] for details on this analytical technique).

In the second region - during the merger - the orbital speed of the binary
gets closer to the speed of light, rendering the PN approximation invalid.
Here the dynamics of the binary is modeled by numerically solving Einstein’s
equations of GR on supercomputers and extrapolating at larger distance to
calculate the resulting waveform (see [34] for an excellent review). Informed
by numerical relativity (NR) [35, 36, 37, 25, 38, 39], the analytical models
are extended to envelop the nonlinear physics of the merger, by using either
the empirical IMRPhenom approach [40, 41], or the EOB reformulation the PN
theory [42, 43, 44]. In the third region, the new black hole formed after
collision settles down by ringing out the extra energy in form of quasi-normal
modes radiation. Here analytical models are valid again and the GWs are
calculated with the black hole perturbation theory [45, 46, 47, 48, 49].

The state-of-art methods mentioned above are continuously improved, driven by
the increasing demands of the GW detection. However, most of them are limited
to quasi-circular orbits (for good reviews see [50, 51]). The LALSuite data
analysis libraries has only one routine for generating eccentric GW templates
based on [52], and current analytical models including eccentricity [52, 53,
54] are under development to increase the accuracy and include all the complex
physics of elliptic binaries [55, 56].

Motivated by the impending detection of GWs from eccentric binaries, we align
our present work with the demand for eccentric waveforms and focus on
accurately implementing fully analytical, easily reproducible GW templates
from eccentric binary black hole (BBH) mergers. We construct a series of new,
fully analytical notebooks within the Wolfram Mathematica software for
calculating the complete GW emitted by compact binaries in eccentric orbits,
by starting with their inspiral, following through the highly nonlinear regime
of the coalescence and ending with the ringdown phase. We describe in detail
how the entire template is built, by bringing together valuable information
from research papers that focus only on specific regions of the orbit and
exhibiting a high level of complexity hidden in scientific jargon, which makes
the problem hard to understand. We open this black box and bring its intricate
content to an approachable level, in a clear way, keeping the balance with the
scientific rigor. Our goal is for the reader to understand, follow and
reproduce our results, even without formal training in GW research.

Another novelty of our implementation resides in the construction of the
hybrid GW strain for the whole evolution of the binary, matching the eccentric
inspiral with the analytic quasi-circular Backward-one-Body (BoB) model [57]
for the merger, by completing them in frequency, amplitude and phase. We find
that for low eccentricity we can connect the two waveforms with no ambiguity
in the time interval, result relevant for the production of accurate GW
waveforms from eccentric binaries. Our research is applicable to the field of
GW astronomy and offers researchers in the field a new tool for generating
eccentric GW, as well as a blueprint for understanding, reproducing and
extending our results, in order to reveal the physics of this kind of sources,
in preparation for their imminent detection.

We start our presentation with a review of the Keplerian problem for eccentric
binaries and its extension to the post-Keplerian parametrization in Section 2.
We clarify the complex physical and mathematical concepts required to solve
analytically the problem of motion when eccentricity is intertwined with
strong gravity, a topic often either overlooked, or hidden in heavy technical
language. Our choice of solving the post-Newtonian Kepler’s equation (based on
[58]) is less used and has the advantage that is purely analytical. We finish
this section by outlining the steps followed to assemble the analytic
gravitational waveform during the inspiral. We continue in Section 3 by
exposing first the difficulty in analytically determining the location of the
transition between the inspiral and merger, then we follow by summarizing the
important formulas employed by the merger model to arrive at the GW strain.

Once the theoretical framework is laid out clearly, it is time to detail our
new implementation and to explain how we built the complete waveform in
Section 4. We start by establishing the time bounds for the binary evolution,
then we implement the inspiral waveform. Lastly, we generate the complete GW
template for the whole evolution of the binary, by matching the GW obtained
for the inspiral with the merger GW in amplitude, frequency and phase, at a
fixed transition time. We employ first small eccentricity, then we push the
boundaries by devising a hybrid waveform for large eccentricity. We conclude
in Section 5 with a summary of our findings, their implications for the field
of gravitational wave astrophysics, and suggest extensions.

We work in geometrized units with $G=c=1$, thus expressing time, space and
mass function the mass of the binary system at infinity
$M=M(\texttt{kg})/M_{\odot}(\texttt{kg})$ where $M_{\odot}$ is the mass of the
sun. This allows an easy rescaling of our results for a given total mass. The
physical units are recovered by multiplying mass with
$M_{\odot}(\texttt{kg})=1.98892\times 10^{30}\texttt{kg}$, space with
$M_{\odot}(\texttt{km})=1.477~{}\texttt{km}$, and time with
$M_{\odot}(\texttt{s})=4.92674\times 10^{-6}\texttt{s}$.

##  2 Theoretical Framework: the Weak-Field Region

###  2.1 The Keplerian Problem

Let us start with an eccentric pair of two black holes reduced to material
points of masses $m_{1}$ and $m_{2}$, separated by a distance $r$, each
revolving in elliptical orbits around a common center of mass. This model is
known as the _two-body problem_. Imposing the conservation of momentum, this
is further simplified to the _one-body problem_ for a single particle of
_reduced mass_ $\mu=m_{1}m_{2}/M$ and position $r$, moving in the external
gravitational field created by the mass $M=m_{1}+m_{2}$ located at the center
of mass. Due to the conservation of angular momentum, the position vector will
be confined to the plane, allowing us to simplify the dynamics to a two-
dimensional problem. Moreover, because $r=r(\cos\theta,\sin\theta)$, we only
need to know how the angle $\theta$ varies in time, to describe the motion of
the point mass $\mu$ [59]. This deceptively simple, now _one-dimensional_
ansatz is also known as the _Keplerina problem_ , but even in Newtonian
gravity it is very complicated to solve. At the heart of this model lies
Kepler’s equation, with its still unknown exact analytical solution, although
scientists tried to find one for centuries. Luckily, we have analytical
methods that approximate this solution with the desired degree of accuracy
[60].

At a a closer look, the elliptical orbit can be described by only four
quantities: the two extreme points that bound the trajectory, $r_{-}$
(periastron, or point of closest approach to the origin) and $r_{+}$
(apoastron or maximum distance from the origin), the semi-major axis
$a=\frac{1}{2}(r_{+}+r_{-})$, and the eccentricity

$e=\frac{r_{+}-r_{-}}{r_{+}+r_{-}}.$ (1)

We align the $x$-coordinate with the major axis of the ellipse and draw a
circle of radius $a$ with origin in the center of the ellipse. Then, we raise
a perpendicular from the x-axis through the orbiting point-mass $\mu$, thus
projecting its motion onto this circle. The motion is now reduced to a
circular orbit on this _reference circle_ described by the angular coordinate
$u$, called _eccentric anomaly_ (see Fig. 1).

Figure 1: Representation of the elliptical orbit as depicted in [61]

In order to calculate the angle $u$, we must use the Kepler’s equation,
deduced from purely geometric considerations as a relationship between this
eccentric anomaly and the eccentricity of the ellipse:

$\ell=u-e\sin u,$ (2)

where $\ell=w(t-t_{0})$ is called the _mean anomaly_ , and $w=2\pi t/P$ is the
_mean motion_ , namely the angular speed of the body $\mu$ averaged over the
period $P$ of the orbit. This quantity is measured with respect to the
periastron, between $\ell=0$ and $\ell=2\pi$ and given by Kepler’s third law,
$w=\sqrt{Ma^{-3}}$ [62]. Once we know this angle, we calculate the separation
with the _shape equation_ $r=a(1-e\cos u)$. The location of each star in the
binary with respect to the center of mass will be recovered with
$|r_{1,2}|=r{m_{1,2}}/M$. Kepler’s equation (2) is transcendental because it
contains a trigonometric function and cannot be computed exactly, although an
exact solution should exist.

Over time, many numerical and analytical approximate solutions have been
proposed, and with the aid of computers we can calculate them with high degree
of accuracy [60]. We will use Bessel’s approach who, using the well known
property that any continuous function can be expanded in Fourier series, wrote
the solution to Kepler’s equation (2) as:

$\displaystyle
u=\ell+\sum_{i=1}^{\infty}a_{i}\sin(i\ell),~{}~{}a_{i}=\frac{2}{i}J_{i}(ie)\sin(i\ell),$
(3)

where $J_{i}(x)$ are the modified Bessel function of the first kind [63]. This
gives an analytical expression for the eccentric anomaly, which is
geometrically related to the orbital phase $\phi$, also known as the _true
anomaly_ , by the equation:

${\phi}=2\tan^{-1}\left(\sqrt{\frac{1+e}{1-e}}\tan\frac{u}{2}\right).$ (4)

Once we know the orbital phase, in Newtonian gravity the problem is solved.

###  2.2 The quasi-Keplerian Problem

Up to this point we have not taken into account the loss of orbital energy due
to the emission of GWs. Let us now include this essential component into our
treatment of the binary’s orbit, while assuming that the orbit decays much
slower compared to the orbital period. This important conjecture, called the
_adiabatic approximation_ , allows us to model the emission of GW with the
_balance equation_ , stating that the loss of orbital energy is balanced by
the flux of GWs for a binary system in a slowly-evolving circular orbit:

$-\frac{dE(t)}{dt}=\dot{E}(t)={\cal F}(t).$ (5)

Here, $E(t)$ is the binding energy of the binary and ${\cal F}(t)$ is the
energy flux emitted in GWs. This equation is valid for most of the inspiral,
as long as the separation $r$ is large and the orbital speed is much smaller
than the speed of light, such that $x_{PN}\ll 1$. In this region, the
gravitational field is weak and the GR equations can be solved in terms of
power series expansions of the parameter $x_{PN}$, starting from the familiar
Newtonian two-body equation of motion. This is the post-Newtonian formalism,
introduced in fact by Einstein [64], and to the reader interested in learning
about it we recommend [65, 32]. This powerful and accurate method provides the
expressions for the binding energy and the emitted flux, enabling us to
integrate the balance equation (5) and to obtain the evolution equation for
the parameter $x_{PN}$

$\frac{dx_{PN}(t)}{dt}=\dot{x}_{PN}(t)=-\frac{{\cal F}(t)}{(dE(t)/dx_{PN})}.$
(6)

The dynamics of the binary becomes much more complicated when we add
eccentricity to the orbit. Due to the ellipticity, the curvature of the orbit
changes, making the angle of the tangent to the orbit vary along the
trajectory, effect described by the _extrinsic curvature_ in GR. This theory
assumes a four-dimensional spacetime, where the orbit depends on how we chose
the three-dimensional spatial slice in this 4D geometry. This means that how
we track the evolution of the eccentricity is not invariant of the coordinate
system, and will depend on how we define it. In mathematical language, the
eccentricity is not _gauge invariant_ and depends on our choice of coordinate
system. To circumvent this problem, we employ again the adiabatic
approximation and assume that the change in eccentricity due to the emission
of GWs occurs on a timescale much longer than the orbital timescale. This
allows us to average the eccentricity over the orbital period and to track its
variation only from one orbit to the next, computing thus the orbit-averaged
version of the evolution equation for the eccentricity. As a note of caution,
this approximation breaks down in the late inspiral of the binary, when the
orbit becomes quasi-circular and the orbital velocity approaches the speed of
light. In this region, starting roughly around $v\approx 0.25$, the emission
of GWs forces the orbit to decay faster and to circularize under the effect of
radiation-reaction [66]. As long as the inspiral does not reach this region,
the orbit decays slowly enough for the motion within one orbit to be
conservative, keeping the adiabatic approximation valid with high degree of
accuracy.

What happens thought when the eccentricity decreases in time due to the
emission of GWs, while the orbital velocity increases? This is the two-body
problem in GR, which does not have an analytic solution. Nevertheless, we can
generalize the solution to Kepler’s equation to include these relativistic
effects, following [58].

$\displaystyle u$ $\displaystyle=\ell+2\sum_{i=1}^{\infty}A_{i}\sin(i\ell),$
(7) $\displaystyle A_{i}$
$\displaystyle=\frac{1}{i}J_{i}(ie_{t})+\sum_{j=1}^{\infty}\alpha_{j}[J_{i+j}(ie_{t})-J_{i-j}(ie_{t})],$
(8)

where $e_{t}$ is the time-dependent eccentricity, $J_{i}$ are again the
modified Bessel functions, and $\alpha_{j}$ are the perturbative, post-
Newtonian corrections, that depend on $e_{t}$, on the parameter $x_{PN}$, as
well as on the chosen coordinate system. This slightly perturbed form of
Kepler’s equation (7) can still accurately describe the orbital evolution in
GR. This formalism is called the _quasi-Keplerian_ parametrization and comes
at a price, namely we must take into account three eccentricity parameters to
include the three important changes in the orbit introduced in GR, namely: (1)
$e_{t}$ for the change in the orbital period due to the emission of GW, (2)
$e_{\theta}$ for the relativistic precession of periastron and (3) $e_{r}$ for
the shrinking of the orbit [67, 68, 69]. These eccentricities depend all on
the energy and angular momentum of the binary and thus on each other, allowing
us to choose in our calculation only one of them, namely $e_{t}$, as is
customary in the literature [11]. Finding analytic solutions to the quasi-
Keplerian equation is, as expected, much more complicated now, but it has been
achieved up to 3PN order in accuracy [58, 70, 71, 72].

###  2.3 The Inspiral Waveform

Up to this point, we described in detail the motion of the binary in the weak-
field region. Armed with this knowledge, we can turn now to our main
objective, namely to calculate the analytical form of the emitted GW. We will
start with the formula for the dimensionless strain,

$h(t)=h_{+}(t)-\mathrm{i}h_{\times}(t)=A(t)e^{-\mathrm{i}\phi_{GW}(t)},$ (9)

where the two strain components ($h_{+}(t)$, $h_{\times}(t)$) describe the two
possible polarizations of the wave, $\phi_{GW}(t)$ is the phase of the GW and
$A(t)$ is its complex amplitude. At the _optimal orientation_ , when the
spherical polar angles of the observer are $\theta=\varphi=0$ [73], the binary
is face-on and right overhead the detector, at some distance $R$ from it. In
this case, the amplitude of GW is maximum, and is given by [53]:

$\displaystyle A(t)$ $\displaystyle=A_{\cal R}(t)+\mathrm{i}A_{\cal I}(t)$
$\displaystyle=-\frac{2\mu}{R}\left(\dot{r}^{2}(t)+r^{2}(t)\dot{\phi}^{2}(t)+\frac{M}{r(t)}+2\mathrm{i}r(t)\dot{r}(t)\dot{\phi}(t)\right).$

GWs are _quadrupolar_ due to the law of conservation of linear and angular
momentum, such that the signal goes through maxima and minima twice during one
orbital cycle, thus: $\phi_{GW}(t)=2\phi(t)$. The strain from eq.(9) is
usually decomposed in spin-weighted ($s=-2$) spherical harmonics modes $(l,m)$
[74],

$h(t)=\sum_{l=2}^{\infty}\sum_{m=-l}^{l}{}_{-2}Y_{lm}h_{lm}(t).$ (11)

The leading mode $(l=2,m=\pm 2)$ typically dominates the sum
($h(t)\approx\sqrt{5/4\pi}h_{22}(t)$).

To calculate the strain in eq.(9) we need to know the equations of motion for
the orbital phase $\phi(t)$ and the orbital separation $r(t)$. In the PN
approximation those quantities are written as powers of the small parameter
$x_{PN}(t)$, as shown below:

$r_{PN}(t)=\frac{M}{x_{PN}(t)}\sum_{j=0}^{N}\rho_{jPN}x_{PN}^{j}(t),$ (12)
$\dot{\phi}_{PN}(t)=\frac{x_{PN}^{3/2}(t)}{M}\sum_{j=0}^{N}\varphi_{jPN}x_{PN}^{j}(t).$
(13)

Here $N$ is the order of the post-Newtonian expansion and the terms
($\rho_{jPN},\varphi_{jPN}$) are numerical expansion coefficients. To
calculate those expressions we need to know how $x_{PN}(t)$ evolves in time,
as given by the balance equation (6). There are several ways of further
expanding this equation in post-Newtonian powers of $x_{PN}(t)$, and we will
choose the _TaylorT4 approximant_ [75], which gives for $x_{PN}(t)$ the
equation [55]:

$\dot{x}_{PN}(t)=\frac{x_{PN}^{5}(t)}{M}\sum_{j=0}^{N}\xi_{jPN}x_{PN}^{j}(t)+\dot{x}_{HT}(t).$
(14)

Note that $\dot{x}_{HT}$(t) are called hereditary terms and enter in the
equation as fractional $j/2$ PN corrections [76]. Those terms describe the
nonlinear interaction of the GWs propagating away from the sources at later
time with the GWs that propagated towards the source at an earlier time
(_memory effects_) and again back-scattered by the spacetime curvature of the
binary (_tail effects_). Those hereditary terms, containing both memory and
tail effects, depend on all the past dynamical history of the source (see [32,
50] for details).

We must add to eq.(14) the following equation of motion for the eccentricity,
because the eccentricity diminishes as the orbit circularizes while the binary
inspirals emitting GWs as it approaches the merger,

$\dot{e}_{PN}(t)=\frac{x_{PN}^{4}(t)}{M}\sum_{j=0}^{N}\epsilon_{jPN}x_{PN}^{j}(t)+\dot{e}_{HT}(t).$
(15)

Next, we find the evolution of the mean anomaly $\ell$ that describes how the
orbital angular velocity changes in the course of an orbital evolution due to
the eccentricity.

$\dot{\ell}_{PN}(t)=w(t)=\frac{x_{PN}^{3/2}(t)}{M}\sum_{j=0}^{N}\lambda_{jPN}x_{PN}^{j/2}(t).$
(16)

After we gather the expansion coefficients
$(\xi_{jPN},\epsilon_{jPN},\lambda_{jPN})$ and calculate how the key
quantities $e_{PN}(t)$, $x_{PN}(t)$ and $l_{PN}(t)$ evolve in time with
eqs.(14), (15) and (16), we must solve the quasi-Keplerian equation (7) to
find the eccentric anomaly $u_{PN}$. The last two steps, before obtaining the
strain in eq.(9), are to integrate equation (13) to find the orbital phase
$\phi(t)$, and to calculate the separation $r(t)$ with equation (12).

##  3 Theoretical Framework: the Strong-Field Region

###  3.1 Location, Location, Location

Let us now attempt to tackle the problem of the binary motion when the orbital
velocity approaches the speed of light and the PN formalism can no longer be
applied.

Going beyond this approximation brings us up against the strong gravitational
field region, which can be described accurately only by Einstein equations of
GR that connect the gravitational field with the geometry of spacetime through
the curvature.

The transition between the weak and strong field starts at the last stable
orbit a particle would have when moving around the central black hole, called
the _innermost stable circular orbit_ (ISCO). This is well defined only for a
small point-like body orbiting a much more massive one, and becomes
increasingly less delimited as the masses of the two bodies are comparable
[77]. We will stay within the one-body problem model considered before, in
which a particle of mass $\mu$ is orbits a central black hole of mass $M$,
only this time we will extend the size of the central body from a material
point to a static Schwarzschild black hole of radius $r_{Sch}$ (the
Schwarzschild radius). In this case, it is analytically proved that ISCO is
located at $r_{ISCO}=6M$ or $3r_{Sch}$ in geometrical units [78]. If the
central black hole is spinning though, ISCO will move closer or further away
from this position, depending on the direction of the particle’s orbit
relative to the black hole’s spin $S_{i}$. The analytical expression for the
dimensionless radius of the last stable prograde orbit is given below [79,
80].

${\tilde{r}_{ISCO}}=\frac{r_{ISCO}}{M}=3+Z_{2}-\sqrt{(3-Z_{1})(3+Z_{1}+2Z_{2})}.$
(17)

Here $Z_{1}=1+(1-\chi_{f}^{2})^{1/3}[(1+\chi_{f})^{1/3}+(1-\chi_{f})^{1/3}]$
and $Z_{2}=\sqrt{3\chi_{f}^{2}+Z_{1}^{2}}$ with $\chi_{f}=S_{f}/M_{f}^{2}$ is
the dimensionless spin of the final black hole.

Beyond ISCO we enter the _plunge_ region, which ends at the _light ring_ (LR)
– the last stable photon orbit. The LR is located at $r_{LR}=1.5r_{Sch}=3M$
and changes if the black hole is spinning, depending on the orbit of the
photon with respect to its spin. Its location, calculated analytically in [79,
81] for a rotating black hole, plays an important role in the matching
procedure, and for the prograde orbit is given by:

$r_{LR}=2M\left[1+\cos{\left(\frac{2}{3}\cos^{-1}(-\chi_{f})\right)}\right].$
(18)

After the LR, the _coalescence_ (or merger) begins, during which the two black
holes come in contact and collide, forming a highly distorted common envelope,
called _apparent horizon_ (AH) [82, 83]. This is a third dimensional surface,
locally defined as the place where the photons directed outward from the
interior region cannot escape to the exterior. The AH coincides with the
boundary of the black hole, namely the _event horizon_ (EH), only for static
black holes $r_{EH}=2M$. This is because the EH is a four dimensional surface
in GR and it’s clear location can be calculated only at the end of all times.
In dynamical spacetimes the EH is approximated with the AH, which is used to
mark the surface of the distorted black hole formed by collision [84]. The
general expression for the radius of the outer event horizon of a rotating
(Kerr) black hole is:

$r_{AH}=r_{+}=M(1+\sqrt{1-\chi_{f}^{2}}).$ (19)

After the merger, we must switch models from the one-body problem, to the
_close limit approximation_ (CLA) [85]. This is an analytical approach that
uses the perturbation theory to analyze how a highly perturbed Schwarzschild
black hole rings down the remaining energy in form of GWs and becomes
quiescent [86, 87]. The amplitude of the emitted GWs increases between the
ISCO and the LR, reaches its peak (or chirps) around the AH, and diminishes in
the ringdown phase, while the frequency increases monotonically.

We separate thus the strong field in three regions: plunge, merger and
ringdown, delimited by three locations: $r_{ISCO}$ (17), $r_{LR}$ (18), and
$r_{AH}$ (18). These locations enter in the analytical calculation of the GWs
emitted in the strong field region and play a key role in determining the
optimal place where we should glue together the weak and strong field
waveforms.

As a side note, it is worth mentioning that several authors [80, 44, 41]
replace the location of the LR with the position of the _minimal energy
circular orbit_ (MECO). This is based on the observation that no constraining
connection exists between the LR, which marks a peak in the effective radial
potential, and the dynamics of the spacetime during the collision, which is
governed by the curvature potential, as is argued in [88]. However, in the
same paper it is pointed out that those two locations are close enough to each
other that the wavelength of the GW is insensitive to such small variation in
the position of the peak. We will consider in this work that the LR is an
accurate approximation for the peak of the curvature potential.

###  3.2 The Merger Waveform

For the merger, we employ the _Backward-one-Body_ (BoB) model, introduced in
[57], that calculates the strain of the GWs using the formula:

$h_{BoB}(t)=-\frac{A_{BoB}(t)}{\Omega^{2}(t)}e^{-i\phi_{BoB}(t)}$ (20)

where the amplitude has the form:

$A(t)_{BoB}=A_{0}\mbox{sech}\left(\frac{t-t_{0}}{\tau_{QNM}}\right)$ (21)

Here $A_{0}$ is an integration constant that enters as a scaling factor,
$t_{0}$ is the time at which the strain of the GW reaches its peak amplitude,
and $\tau_{QNM}$ is the damping time necessary for the final black hole to
settle down after the collision. During this time, the highly perturbed
remnant object releases its energy in GWs as characteristic exponentially
damped _quasi-normal_ modes (QNM). We can find the damping time using the
definition [89]:

$\tau_{QNM}=2\frac{Q_{QNM}}{\Omega_{QNM}}.$ (22)

where the quality factor $Q_{QNM}$ is a dimensionless quantity measuring the
number of oscillations observed before a mode gets attenuated by a factor of
$e^{2\pi}$ [90], and $\Omega_{QNM}$ is the QNM frequency. Whatever initial
data we provide inside the $r<3M$ region where the peak of the potential
barrier is located, this barrier will filter it, and an outside observer will
detect only the QNM ringing. We employ the following polynomial fit for the
pair ($\Omega_{QNM},Q_{QNM}$):

$M_{f}\Omega_{QNM}=f_{1}+f_{2}(1-\chi_{f})^{f_{3}},$ (23)
$Q_{QNM}=q_{1}+q_{2}(1-\chi_{f})^{q_{3}}$ (24)

with the set of coefficients $f_{1,2,3}$ and $q_{1,2,3}$ given in [48].

The phase in eq.(20) is obtained by integrating the frequency $\Omega(t)$:

$\Phi_{BoB}(t)=\int\Omega_{BoB}(t)dt,$ (25)

which is given in this model by the equation:

$\Omega_{BoB}(t)=\left(\Omega_{i}^{4}+\kappa\left[\tanh\left(\frac{t-t_{0}}{\tau}\right)-\tanh\left(\frac{t_{i}-t_{0}}{\tau}\right)\right]\right)^{1/4}$
(26)

This equation makes use of the parameter $k$ given by the expression:

$\kappa=\left[\frac{\Omega_{QNM}^{4}-\Omega_{i}^{4}}{1-\tanh\left(\frac{t_{i}-t_{0}}{\tau}\right)}\right]$
(27)

The initial frequency $\Omega_{i}$ corresponds to the frequency at the end of
the inspiral, around the LR, where the PN approximation breaks down. The most
important quantity is the initial time $t_{i}$, which is analytically obtained
by replacing eq.(27) in eq.(26):

$t_{i}=t_{0}-\frac{\tau_{QNM}}{2}\ln\left(\frac{\Omega_{QNM}^{4}-\Omega_{i}^{4}}{2\tau\Omega_{i}^{3}\dot{\Omega}_{i}}-1\right).$
(28)

This equation is essential, being the key ingredient that determines the exact
time when the merger waveform is added to the inspiral when constructing the
total strain.

##  4 Results: Building the Compete Waveform

###  4.1 It’s about time

Before we turn all this formalism into concrete waveforms, let us take a
better look at the independent variable dictating the evolution of the whole
system, namely the time coordinate, which is taken to coincide with the proper
time of distant observers.

Let’s start with the expression given in [7] for the _coalescence time_ until
a binary in quasi-circular orbit reaches the merger, from a known initial
separation or semi-major axis $r_{0}$, due to the emission of GW. In
geometrical units, this formula is

$T_{c}=\frac{5}{256}\frac{r^{4}_{0}}{m_{1}m_{2}(m_{1}+m_{2})}=\frac{5}{256}\frac{r^{4}_{0}(1+q)^{2}}{M^{3}q}.$
(29)

where $q=m_{1}/m_{2}$ is the mass ratio of the binary. The addition of a small
initial eccentricity $e_{0}$ is accounted for by replacing $r_{0}$ in eq.(29)
with

$a_{0}=r_{0}\frac{8.70127(1-e_{0}^{2})}{(304+121e_{0}^{2})^{870/229}}.$ (30)

Recently, improved estimations of this time, including $1^{st}$ order PN
perturbations, were given in [91, 92]. The choice of the initial separation
$a_{0}$ is limited by the low frequency cutoff of Earth’s seismic activity
$f_{low}$, which dictates the threshold value of GWs detection band for aLIGO.
For a cutoff frequency $f_{low}=10~{}\texttt{Hz}$, we obtain from Kepler’s
third law $a_{0}=({M}/{\pi^{2}f_{low}^{2}})^{1/3}$. For an equal mass binary
of total mass $M=50M_{\odot}$, this amounts to a separation of about
$1278M_{\odot}=1888~{}\texttt{km}$, and gives a coalescence time
$T_{c}=8.22~{}\texttt{s}$ of the GW signal in the aLIGO band. The coalescence
time decreases for eccentric binaries, thus if we assume for example that the
binary enters the aLIGO band with $e_{0}=0.25$, it will take
$T_{c}=6.11~{}\texttt{s}$ to collide, which will be enough time for the
eccentricity to be detected. The LISA band cutoff frequency however is much
lower, $f_{low,LISA}=10^{-4}~{}\texttt{Hz}$, which, for a $50M_{\odot}$ equal-
mass binary gives a separation of $4\times 10^{6}~{}\texttt{km}$, and
consequently a collision time of about five millions years. Thus systems with
large eccentricities are expected to be regularly detected by space-borne
instruments. Conversely, GW signals lasting only a few seconds in the LISA
band will come from billions solar masses binaries.

The issue we must tackle next is to find an appropriate estimate for the time
when the inspiral model breaks down and must be replaced with the merger
model. As we mentioned before, the transition between the weak and strong
regime happens around the last stable orbit (ISCO), located around
$r_{ISCO}=6M$. The problem with the start of the merger phase at ISCO is, as
we explained previously in Sec.3.1, that ISCO is not well defined for
comparable mass binaries. Fortunately, we can still define an approximate ISCO
even in this case, closer to the light-ring, given in eq.(18), which in the
non-spinning case ($\chi_{f}=0$) is $r_{LR}=4M$ [43]. Using eq.(29) for
$r_{0}=r_{LR}=4M$, we get $t=20M$, which corresponds to the initial time
$t_{i}$ for the strong region. Here, in the $0^{th}$ PN approximation, the
orbital frequency will be $M\omega_{LR}=0.125$, which represents a lower limit
for the orbital frequency of the inspiral model predicted by the PN theory
around $t_{LR}$. Another constraint, offered by eq.(28), provides a third
independent check for the relationship between this time and the frequency at
which the merger waveform joins smoothly with the inspiral. The reader is
warned that some arbitrariness will be still inherent due to the assumptions
we will make, because there is no unique way of defining this transition time
unambiguously.

Lastly, we will choose the time at which the amplitude of the fundamental mode
has a maximum, as the origin of the time axis $t_{0}=0$.

###  4.2 Eccentric Slowness

Our procedure for the implementation of the inspiral waveform follows the
outline given in Sec.2.3. First, we implement all the expressions for the
coefficients that appear in equations (12), (13), (14), (15) and (16). Those
coefficients are currently known up to the $4^{th}$ leading PN order, with up
to 6PN quasi-circular correction terms and hereditary corrections [93, 94, 95,
33]. Impressive efforts are continuously made to calculate new terms in the PN
expansion [96, 97, 98, 99, 100]. Following [76], we include up to 6PN self-
force and hereditary correction terms, to increase the accuracy in modeling
the region near the merger. We use NDSolve [101] within the Wolfram
Mathematica symbolic computation program [102] to solve the coupled equations
of motion (14) and (15), with the $\xi_{jPN}$ terms as in [76] and the
eccentricity coefficients $\epsilon_{j_{PN}}$ from [73, 103]. We start from
the aLIGO cutoff frequency $f_{low}$, work with $20$ precision digits and pick
the collisions time $T_{s}$ when the equations become stiff and fail to yield
a solution, then mark $t_{LR}=T_{s}-20M$ as the end of the weak field region.

Let’s look now at the predicted evolution of the eccentricity, for an equal
mass binary of $50M_{\odot}$ total mass and three values for the eccentricity:
$e_{0}=0.5,0.25,0.1$ at $f_{low}$. We choose the highest eccentricity to be
$e_{0}=0.5$ because is close to the average eccentricity for wide binaries
[104]. Fig.2 shows that collisions of highly eccentric binaries happen faster,
thus we expect a higher number of GW detections to come from initially highly
eccentric systems.

Figure 2: Time evolution of the eccentricity, for an equal mass binary of
$50M_{\odot}$ (the time is measured in units of mass, $1M=4.92674\times
10^{-6}s$).

Fig. 3, shows the effect of the mass ratio on the eccentricity evolution,
starting from $e_{0}=0.25$ eccentricity at cutoff frequency, for a
$50M_{\odot}$ mass binary with three mass ratios $q=4,2,1$. We see that small-
mass ratio binaries take less time to merge than high-mass ratio ones, which
has the expected consequence that more GW signals will come from equal or
nearly equal mass binaries.

Figure 3: Time evolution of the eccentricity, for a mass binary of
$50M_{\odot}$ and three mass ratios: $q=1,2,4$.

As expected, the evolution of the eccentricity depends also on the method
used. Fig. 4 shows that taking into account the hereditary terms is important,
because they add the fractional $n/2$ powers and cancel out terms depending on
the chosen system of coordinates, thus ensuring the _gauge-invariance_ of the
equations [76].

Figure 4: Time evolution of the eccentricity, for a mass binary of
$50M_{\odot}$, mass ratio $q=2$ and eccentricity $e_{0}=0.25$, with three PN
models.

In Table 1 we give the values for the eccentricity, slowness, and orbital
frequency at the LR for an equal-mass binary with total mass $50$M, starting
from three different eccentricities, as well as a comparison between the
numerical value of the collision time obtained when the equations become stiff
$T_{s}$ and the collision time $T_{c}$ predicted by eq.(29). The slowness
reaches a constant value at the LR ($x_{LR}\approx 0.178$), while the
eccentricity decreases dramatically, more than $15$ times near LR, which
proves that by the beginning of the strong field region the orbit of the
system is circularized, even for large initial eccentricity, which can happen
for capture binaries in dense clusters [105]. It is reassuring that the
eccentricity becomes close enough to zero by the time we reach the strong-
field region, because the merger model is built for quasi-circular binaries.

Table 1: Comparison between the computed time to LR and
the estimated collision time for three initial eccentricities at $f_{low}$,
$x_{PN}$ and the angular frequency at the LR for a $50M_{\odot}$ binary.

  * $e_{0}$ | $e_{LR}$ | $x_{LR}$ | $M\Omega_{LR}$ | $T_{LR}(s)$ | $T_{c}(s)$
---|---|---|---|---|---
$0.10$ | $5.952\times 10^{-3}$ | $0.1782$ | $7.525\times 10^{-2}$ | $7.746$ | $7.842$
$0.25$ | $1.671\times 10^{-2}$ | $0.1782$ | $7.526\times 10^{-2}$ | $6.074$ | $6.113$
$0.50$ | $5.779\times 10^{-2}$ | $0.1776$ | $7.483\times 10^{-2}$ | $2.010$ | $2.251$

###  4.3 All kinds of anomalies

In a previous work [106] we implemented an analytical calculation of GW from
quasi-circular binaries up to the $3.5PN$ order, first by solving only
eq.(14), then by relying on Kepler’s third law to obtain the phase as:
$\dot{\phi}_{PN}(t)=\omega_{PN}(t)=M^{-1}x_{PN}^{3/2}(t)$. The addition of the
eccentricity complicates things by introducing the angular anomaly which
requires solving the transcendental Kepler’s equation.

The orbital angular velocity is now represented by the mean motion $w$, which
changes in the course of an orbital evolution due to the presence of the
eccentricity and is described by eq.(16). To find the evolution of the mean
anomaly $l(t)$ we integrate this equation with the $\lambda_{jPN}$
coefficients from [73]. With this quantity in hand, we proceed to calculating
the eccentric anomaly given in eq.(7), with the $3$PN coefficients
$\alpha_{j}$ entering in the term $A_{i}$ taken from [58] for the Arnowitt-
Deser-Misner (ADM) coordinates. The evaluation of the Bessel functions is
computationally expensive and in order to speed up the calculation we devised
a method detailed in [107]. In here, we increase the truncation of the sum in
the Bessel functions to $10$, tabulate $u(t)$, then apply a high order
interpolating polynomial to keep the error around $10^{-20}$ precision.

With the hardest part of the implementation out of our way, the next two steps
are to calculate the orbital separation $r(t)$ given by eq.(12) and the
derivative of the phase as in eq.(13), with the coefficients $\rho_{jPN}$ and
$\phi_{jPN}$ from [73]. In order to keep the desired precision and speed up
the calculation, before we evaluate the time derivative of the separation, we
use the same procedure of fitting a high-order polynomial to $r(t)$. Lastly,
we solve for the phase $\phi(t)$ and build the strain with eq.(9). We use the
following formulas for the polarizations of the strain for optimal orientation
of the observer [106]:

$\displaystyle h_{+}(t)=-2\frac{M\eta}{R}\Bigg{[}\Bigg{.}$
$\displaystyle\left(-\dot{r}(t)+r^{2}(t)\dot{\phi}^{2}(t)+\frac{M}{r(t)}\right)\cos
2\phi(t)$ $\displaystyle+2r(t)\dot{r}(t)\dot{\phi}(t)\sin
2\phi(t)\Bigg{.}\Bigg{]}$ $\displaystyle
h_{\times}(t)=-2\frac{M\eta}{R}\Bigg{[}\Bigg{.}$
$\displaystyle\left(-\dot{r}(t)+r^{2}(t)\dot{\phi}^{2}(t)+\frac{M}{r(t)}\right)\sin
2\phi(t)$ $\displaystyle-2r(t)\dot{r}(t)\dot{\phi}(t)\cos
2\phi(t)\Bigg{.}\Bigg{]}$

Those expressions can be extended to take into account the orientation of the
binary with respect to the detector by performing three consecutive rotations
from the frame of the orbit $\\{x,y,z\\}$ onto the frame of the observer
$\\{X,Y,Z\\}=P\\{x,y,z\\}$, where $P$ is the rotation matrix written function
of the three Euler angles, namely the longitude of the ascending node $\psi$,
the inclination $\iota$ and the argument of the pericenter $\nu$, as described
in [62]. In this case, the true anomaly must be corrected to:
$\phi\rightarrow\Phi=\phi-\nu$. We assume an optimal orientation of the
detector and leave the inclusion of arbitrary orientation for a future work.

We build our implementation hierarchically, constructing first the 2PN
corrections, then adding the 3PN terms, and lastly the hereditary term
$\dot{x}_{HT}$, while carefully testing our results against the zero
eccentricity limit using the first GW detected signal GW150914 as a sanity
check (as detailed in [107]). The addition of the higher order hereditary
corrections is necessary because, besides adding accuracy, it extends the
validity of the inspiral model closer to the merger, increasing the overlap
region with the merger model.

We show in Figures 5 and 6 the amplitude and the argument of the strain for
optimal orientation of the observer. It is clear again, from Fig.5, that the
eccentricity plays a vital role in determining the strength and duration of
the signal, because highly eccentric system are more energetic but much
shorter. In Fig.6 we plot the argument of the strain, which in absolute value
is equal to twice the orbital phase $\phi$. We see that both the amplitude and
the argument (phase) exhibit oscillations that diminish with the decrease of
the eccentricity, until the orbit becomes quasi-circular near the merger.

Figure 5: Time evolution for the amplitude of the GW, with the time measured
in units of mass. Figure 6: Time evolution of the unwrapped argument of the
strain, with time measured in units of mass.

###  4.4 Prepare to Merge

Before we proceed to model the merger part of the GW strain, we must know
first the mass, spin and ringdown frequency of the final black hole, required
by the BoB model as initial data. For the building of the inspiral strain we
have not taken this spin into account because we treated the black holes as
particles. Once the separation becomes comparable to the LR radius, we cannot
ignore the spin any longer. If the spins are equal and antiparallel, or the
black holes don’t spin, then the final spin is $\chi_{0}=0.686$. In this case
eq.(18) shows that the location of the LR will move to about
$r_{LR,\chi_{0}}=2.04M$, and the position of the apparent horizon, given by
eq.(19), will be around $r_{AH,\chi_{0}}=1.73M$.

With the final spin known, we proceed next to calculate the mass loss to GW
energy, which subtracted from the total mass gives the mass of the final black
hole. A simple estimate for the final mass is given by the Christodoulou
formula [108]

$M_{f}=\sqrt{M_{irr}^{2}+\frac{S_{f}^{2}}{4M_{irr}^{2}}},~{}\textrm{where}~{}M_{irr}=\sqrt{\frac{A_{AH}}{16\pi}}=\frac{r_{AH}}{2}.$
(33)

The last essential ingredients we must add to the initial data for the merger
model is the knowledge of fundamental QNM frequency $\Omega_{QNM}$ of the
remnant black hole and the quality factor $Q_{QNM}$, which is necessary to
calculate the damping time describing the decay in the amplitude of the GW.
For this we use eq.(23), (24) with the set of coefficients for the dominant
mode taken from [109]:

$f_{1}=1,~{}f_{2}=-0.63,~{}f_{3}=0.3,~{}g_{1}=0,~{}g_{2}=2,~{}g_{3}=-0.45.$
(34)

The damping time of the fundamental mode is given by eq.(22).

We calculate the final mass $M_{f}$, the dominant $(l=m=2)$ mode of the
resonant frequency $\omega_{QNM}$, the corresponding quality factor $Q_{QNM}$
and the damping time $\tau$ for three non-spinning binary systems with mass
rations $q=1,2,4$ and give the values in Table 2. We add to this table the PN
values for $\Omega_{i}$ and its time derivative $\dot{\Omega}_{i}$, calculated
in the inspiral model at the light ring.

Table 2: Initial data required for the BoB model.

  * $q$ | $\chi_{f}$ | $M_{f}$ | $\omega_{QNM}$ | $Q_{QNM}$ | $\tau$ | $\Omega_{i}$ | $\dot{\Omega}_{i}$
---|---|---|---|---|---|---|---
$1$ | $0.6865$ | $0.9516$ | $0.5698$ | $3.301$ | $11.586$ | $7.529\times 10^{-2}$ | $8.577\times 10^{-4}$
$2$ | $0.6231$ | $0.9612$ | $0.5385$ | $3.056$ | $11.349$ | $7.414\times 10^{-2}$ | $7.399\times 10^{-4}$
$4$ | $0.4637$ | $0.9779$ | $0.4782$ | $2.642$ | $11.049$ | $7.169\times 10^{-2}$ | $5.046\times 10^{-4}$

The BoB model (as emphasized in [57]) does not depend directly on any other
numerically fitted coefficients. We start building it now, by calculating the
angular frequency of the fundamental mode using eq.(26), with $\kappa$ given
by eq.(27) and $t_{i}$ from eq.(28), then we obtain the dominant mode by
multiplying it with $m=2$. Once the frequency is known, we calculate the
amplitude, given by eq.(21). We proceed next to calculate the phase, given by
eq.(25). Lastly, we implement the strain, given by (20). We remind the reader
that we assumed the black holes entering the merger do not spin, or their
effective spin is zero, and indeed this is a good approximation for most of
the binaries detected in the LIGO and Virgo O1–O3 observing runs [4].

###  4.5 Fasten It Tight

It is due time now to attach the GW template for the inspiral to the one
obtained with the BoB mode and to construct a completely analytical GW
template that encompasses the whole evolution of the BBH collision, from the
time it enters the detection band untill the quiescence of the final BH.

We use an equal-mass binary of normalized mass $m_{1}+m_{2}=1M$ and start the
inspiral at a separation $r_{0}=10M$, but those assumptions can be easily
modified to fit any separation, and rescaled to a desired total mass of the
BBH system, as well as mass ratio up to $q=10$. For the initial eccentricity
we will consider two values, $e_{i,1}=10^{-3}$ corresponding to the quasi-
circular case, and an eccentric close orbit of $e_{i,2}=0.15$. With this
initial data we run the PN portion of the binary evolution. The quasi-circular
case evolves up to $T_{s,1}=530M$, its eccentricity drops to
$e_{f,1}=1.7\times 10^{-4}$ and its orbital frequency reaches
$\Omega_{f,1}=0.13407$. We take the initial frequency for the BoB model at
$t_{LR,1}=T_{s,1}-20M$ corresponding to the location of the LR, where the PN
approximation breaks. The value of this frequency, taken directly from the PN
data, is $\Omega_{i,1}=7.4085\times 10^{-2}$. The large eccentricity case has
a shorter span, only up to $T_{s,2}=460M$, its eccentricity drops to
$e_{f,2}=1.988\times 10^{-2}$ and its orbital frequency raises to
$\Omega_{f,2}=0.19466$. Even for this case, the initial frequency for the BoB
model at $t_{LR,2}=T_{f,2}-20M$ is $\Omega_{i,2}=7.5837\times 10^{-2}$, which
is close to the values found in Table 1. Note that the orbital frequency at
the LR stays $\approx 0.075$ with a mean deviation of $\approx 10^{-4}$,
marking the transition from the weak (far field) to the strong (near field)
zone where the reduced wavelength of the gravitational wave
$\bar{\lambda}_{GW}=\Omega^{-1}$ becomes comparable to the binary separation
$r$.

Save for the slight difference in the initial frequency, we start the BoB
model with the same values for the ringdown frequency and quality factor. We
perform the stitching in frequency at the initial time calculated in eq.(28)
$t_{i}$ for the BoB model.

We create the hybrid frequency with the formula:

$\omega_{hyb}(t)=(1-\sigma(t))\omega_{PN}(t-ti)+\sigma(t)\omega_{BoB}(t).$
(35)

where $\sigma(t)$ is the piecewise function:

$\sigma(t)=\left\\{\begin{array}[]{lr}0&t\leq t_{i}\\\
1&t>t_{i}\end{array}\right.$ (36)

We obtain an excellent overlap at and around $t_{i}$, as shown in Fig.7,
without introducing any other time shift. Even the high-eccentricity case, as
shown in Fig.8, exhibits this excellent behavior at the chosen matching time.
We drastically reduced the ambiguity in time for the matching in frequency,
which is a remarkable result.

Figure 7: Hybrid frequency, by matching the PN at $t_{LR,1}$ with BoB at
$t_{i}$, for an equal-mass binary merger with initial eccentricity
$e_{i}=10^{-3}$ at $r=10M$. Figure 8: Hybrid frequency, by matching the PN at
$t_{LR,2}$ with BoB at $t_{i}$, for an equal-mass binary merger with initial
eccentricity $e_{i}=0.15$ at $r=10M$.

Before stitching together the strain, we normalize the BoB amplitude and
rescale the amplitude of the inspiral with this normalized amplitude divided
by the PN amplitude at the light ring $t_{LR}$.

We form the hybrid strain with the same technique, this time translating also
the BoB strain with the peak time of increase in frequency $t_{f}$ and
correcting it for sign, which indicates a phase difference of $\pi/2$. We note
that this does not introduce a shift ambiguity because the two models might
not to use the same sign convention. The formula for the normalized hybrid
strain is:

$\bar{h}_{hyb}(t)=(1-\sigma(t))\bar{h}_{PN}(t-ti)-\sigma(t)\bar{h}_{BoB}(t+t_{p}).$
(37)

For the low-eccentricity case this provides again an excellent overlap at and
around $t_{i}$, as seen in Fig.9.

Figure 9: Hybrid strain, by matching the PN at $t_{LR,1}$ with BoB at $t_{i}$,
translated with $t_{p}$, for an equal-mass binary merger with initial
eccentricity $e_{i}=10^{-3}$ at $r=10M$.

When matching the high-eccentricity inspiral with the BoB merger strain we
encounter a slight mismatch, both in amplitude and time of the matching. This
is most likely due to the fact that the orbit is not fully circularized by the
time the binary enters the strong-field zone, and the BoB model does not take
into account the eccentricity. We mitigate residual effects or the
eccentricity by matching the two strains at the closest peak in their
amplitude near $t_{i}$ and rescaling slightly the inspiral strain to match the
amplitude of the BoB at that time. The result is shown in Fig.10. This indeed
introduces an expected ambiguity in the matching interval for high-
eccentricities, which can be mitigated in further studies.

Figure 10: Hybrid strain, by matching the PN at $t_{LR,2}$ with BoB at
$t_{i}$, translated with $t_{p}/3$, for an equal-mass binary merger with
initial eccentricity $e_{i}=0.15$ at $r=10M$.

##  5 Conclusions

In this work we implemented a purely analytical framework for obtaining
complete GW templates from eccentric binary sources and produced a set of new,
freely available Wolfram Mathematica notebooks. First, we gathered and
assembled into a coherent way all the pieces of the puzzle required to ensure
the accurate implementation. We started with an exposition of the essential
theoretical framework, by introducing the specific terminology and by giving
clear explanations to prevent confusing readers with different scientific
background. We carefully chose, from the rich scientific literature, a fully
analytical procedure to build the inspiral GW that contains high-order PN
corrections for enhanced accuracy. For the eccentricity and slowness we employ
instantaneous and hereditary coefficients with energy corrections up to $6$PN.
We implemented a purely analytical expression for the solution to Kepler’s
equation that gives the eccentric anomaly, accurate up to $3$PN order. We
calculated the orbital separation and the phase, keeping the errors below
machine precision, and used high order interpolating polynomials when
necessary to speed up the calculation. We assembled the strain and displayed
its dependence on eccentricity and mass ratio, showing that we expect a higher
number of GW detections from eccentric comparable mass ratios binaries, which
emphasizes the relevance of modeling such systems. For the merger, we
implemented the purely analytical BoB model for the merger, which relies on
the quasi-circular assumption. At last, we glued together the eccentric
inspiral and the quasi-circular merger waveforms by completing them in
frequency, amplitude and phase with a piecewise function, building the hybrid
GW strain for the whole evolution of the binary. Remarkably, for low
eccentricity, the match between the BoB and PN approximations are on-point,
with coincidence in the overlap, which indicates no ambiguity in the time
interval, a great improvement from the usual matching techniques. For high-
eccentricity we must make a slight adjustment, both in time and amplitude, to
compensate for the implicit quasi-circular assumption built into the BoB
approach, but the matching in frequency and phase are again unambiguous.

The new open-source tool for the calculation of analytic GW templates for
eccentric binaries, together with the thorough and streamlined documentation
of our steps given in this paper, offers researchers in the field of
gravitational waves a straightforward path to understand, reproduce, use and
extend our implementation. Future developments will incorporate higher PN
orders, will include spin, precession and kicks, and might augment the BoB
model to account for eccentricity close to the merger. Our notebooks are
available at github.com/mbabiuc/MathScripts.

The authors wish to acknowledge assistance from the Physics Department at
Marshall University, and financial support from the National Science
Foundation through the EPSCoR Grant OIA-1458952 to the state of West Virginia
“Waves of the Future”.

## References

## References

  * [1] Celoria M, Oliveri R, Sesana A and Mapelli M 2018 Lecture notes on black hole binary astrophysics (Preprint http://arxiv.org/abs/1807.11489)
  * [2] Abbott B P et al (LIGO Scientific, Virgo) 2016 Phys. Rev. Lett. 116 061102 (Preprint http://arxiv.org/abs/1602.03837)
  * [3] Schutz B F 1999 Classical and Quantum Gravity 16 A131–A156 URL https://doi.org/10.1088/0264-9381/16/12a/307
  * [4] Abbott R et al (LIGO Scientific, Virgo) 2021 SoftwareX 13 100658 (Preprint http://arxiv.org/abs/1912.11716)
  * [5] Brown D A, Vahi K, Taufer M, Welch V and Deelman E 2021 Comput. Sci. Eng. 23 73–82 (Preprint http://arxiv.org/abs/2010.07244)
  * [6] Nitz A H, Lenon A and Brown D A 2019 Astrophys. J. 890 1 (Preprint http://arxiv.org/abs/1912.05464)
  * [7] Peters P C 1964 Physical Review 136 1224–1232
  * [8] Petrovich C and Antonini F 2017 Astrophys. J. 846 146 (Preprint http://arxiv.org/abs/1705.05848)
  * [9] Antonini F and Rasio F A 2016 Astrophys. J. 831 187 (Preprint http://arxiv.org/abs/1606.04889)
  * [10] Antonini F, Toonen S and Hamers A S 2017 Astrophys. J. 841 77 (Preprint http://arxiv.org/abs/1703.06614)
  * [11] Loutrel N, Liebersbach S, Yunes N and Cornish N 2019 Classical and Quantum Gravity 36 025004 (Preprint http://arxiv.org/abs/1810.03521)
  * [12] Kyutoku K and Seto N 2014 Mon. Not. Roy. Astron. Soc. 441 1934–1942 (Preprint http://arxiv.org/abs/1312.2953)
  * [13] Hall E D et al 2021 Phys. Rev. D 103 122004 (Preprint http://arxiv.org/abs/2012.03608)
  * [14] Maggiore M et al 2020 JCAP 03 050 (Preprint http://arxiv.org/abs/1912.02622)
  * [15] Baker J et al 2019 Space Based Gravitational Wave Astronomy Beyond LISA (Preprint http://arxiv.org/abs/1907.11305)
  * [16] Barausse E et al 2020 Gen. Rel. Grav. 52 81 (Preprint http://arxiv.org/abs/2001.09793)
  * [17] Isoyama S, Nakano H and Nakamura T 2018 PTEP 2018 073E01 (Preprint http://arxiv.org/abs/1802.06977)
  * [18] Mei J et al (TianQin) 2021 PTEP 2021 05A107 (Preprint http://arxiv.org/abs/2008.10332)
  * [19] Lower M E, Thrane E, Lasky P D and Smith R 2018 Phys. Rev. D 98 083028 (Preprint http://arxiv.org/abs/1806.05350)
  * [20] Gondán L, Kocsis B, Raffai P and Frei Z 2018 Astrophys. J. 855 34 (Preprint http://arxiv.org/abs/1705.10781)
  * [21] Breivik K, Rodriguez C L, Larson S L, Kalogera V and Rasio F A 2016 Astrophys. J. Lett. 830 L18 (Preprint http://arxiv.org/abs/1606.09558)
  * [22] Buonanno A and Damour T 1999 Phys. Rev. D 59 084006 (Preprint http://arxiv.org/abs/gr-qc/9811091)
  * [23] Bohé A et al 2017 Phys. Rev. D 95 044028 (Preprint http://arxiv.org/abs/1611.03703)
  * [24] Nagar A et al 2018 Phys. Rev. D 98 104052 (Preprint http://arxiv.org/abs/1806.01772)
  * [25] Santamaria L et al 2010 Phys. Rev. D 82 064016 (Preprint http://arxiv.org/abs/1005.3306)
  * [26] Estellés H, Ramos-Buades A, Husa S, García-Quirós C, Colleoni M, Haegel L and Jaume R 2021 Phys. Rev. D 103 124060 (Preprint http://arxiv.org/abs/2004.08302)
  * [27] García-Quirós C, Colleoni M, Husa S, Estellés H, Pratten G, Ramos-Buades A, Mateu-Lucena M and Jaume R 2020 Phys. Rev. D 102 064002 (Preprint http://arxiv.org/abs/2001.10914)
  * [28] Pratten G, Husa S, Garcia-Quiros C, Colleoni M, Ramos-Buades A, Estelles H and Jaume R 2020 Phys. Rev. D 102 064001 (Preprint http://arxiv.org/abs/2001.11412)
  * [29] Asada H and Futamase T 1997 Prog. Theor. Phys. Suppl. 128 123–181 (Preprint http://arxiv.org/abs/gr-qc/9806108)
  * [30] Blanchet L 2011 Fundam. Theor. Phys. 162 125–166 (Preprint http://arxiv.org/abs/0907.3596)
  * [31] Will C M 2011 Proc. Nat. Acad. Sci. 108 5938 (Preprint http://arxiv.org/abs/1102.5192)
  * [32] Blanchet L 2014 Living Rev. Rel. 17 2 (Preprint http://arxiv.org/abs/1310.1528)
  * [33] Blanchet L 2019 Comptes Rendus Physique 20 507–520 (Preprint http://arxiv.org/abs/1902.09801)
  * [34] Baumgarte T W and Shapiro S L 2010 Numerical Relativity: Solving Einstein’s Equations on the Computer (Cambridge University Press)
  * [35] Buonanno A, Cook G B and Pretorius F 2007 Phys. Rev. D 75 124018 (Preprint http://arxiv.org/abs/gr-qc/0610122)
  * [36] Hannam M, Husa S, Sperhake U, Bruegmann B and Gonzalez J A 2008 Phys. Rev. D 77 044020 (Preprint http://arxiv.org/abs/0706.1305)
  * [37] Boyle M, Brown D A, Kidder L E, Mroue A H, Pfeiffer H P, Scheel M A, Cook G B and Teukolsky S A 2007 Phys. Rev. D 76 124038 (Preprint http://arxiv.org/abs/0710.0158)
  * [38] Szilágyi B, Blackman J, Buonanno A, Taracchini A, Pfeiffer H P, Scheel M A, Chu T, Kidder L E and Pan Y 2015 Phys. Rev. Lett. 115 031102 (Preprint http://arxiv.org/abs/1502.04953)
  * [39] Hinder I, Ossokine S, Pfeiffer H P and Buonanno A 2019 Phys. Rev. D 99 061501 (Preprint http://arxiv.org/abs/1810.10585)
  * [40] Ajith P et al 2011 Phys. Rev. Lett. 106 241101 (Preprint http://arxiv.org/abs/0909.2867)
  * [41] Jiménez-Forteza X, Keitel D, Husa S, Hannam M, Khan S and Pürrer M 2017 Phys. Rev. D 95 064024 (Preprint http://arxiv.org/abs/1611.00332)
  * [42] Pan Y, Buonanno A, Baker J G, Centrella J, Kelly B J, McWilliams S T, Pretorius F and van Meter J R 2008 Phys. Rev. D 77 024014 (Preprint http://arxiv.org/abs/0704.1964)
  * [43] Buonanno A, Pan Y, Baker J G, Centrella J, Kelly B J, McWilliams S T and van Meter J R 2007 Phys. Rev. D 76 104049 (Preprint http://arxiv.org/abs/0706.3732)
  * [44] Babak S, Taracchini A and Buonanno A 2017 Phys. Rev. D 95 024010 (Preprint http://arxiv.org/abs/1607.05661)
  * [45] Mino Y, Sasaki M, Shibata M, Tagoshi H and Tanaka T 1997 Prog. Theor. Phys. Suppl. 128 1–121 (Preprint http://arxiv.org/abs/gr-qc/9712057)
  * [46] Kokkotas K D and Schmidt B G 1999 Living Rev. Rel. 2 2 (Preprint http://arxiv.org/abs/gr-qc/9909058)
  * [47] Sasaki M and Tagoshi H 2003 Living Rev. Rel. 6 6 (Preprint http://arxiv.org/abs/gr-qc/0306120)
  * [48] Berti E, Cardoso V and Starinets A O 2009 Class. Quant. Grav. 26 163001 (Preprint http://arxiv.org/abs/0905.2975)
  * [49] Harms E, Bernuzzi S, Nagar A and Zenginoglu A 2014 Class. Quant. Grav. 31 245004 (Preprint http://arxiv.org/abs/1406.5983)
  * [50] Le Tiec A 2014 Int. J. Mod. Phys. D 23 1430022 (Preprint http://arxiv.org/abs/1408.5505)
  * [51] Isoyama S, Sturani R and Nakano H 2020 Post-Newtonian templates for gravitational waves from compact binary inspirals (Preprint http://arxiv.org/abs/2012.01350)
  * [52] Tanay S, Haney M and Gopakumar A 2016 Phys. Rev. D 93 064031 (Preprint http://arxiv.org/abs/1602.03081)
  * [53] Hinder I, Kidder L E and Pfeiffer H P 2018 Phys. Rev. D 98 044015 (Preprint http://arxiv.org/abs/1709.02007)
  * [54] Susobhanan A, Gopakumar A, Hobbs G and Taylor S R 2020 Phys. Rev. D 101 043022 (Preprint http://arxiv.org/abs/2002.03285)
  * [55] Huerta E A et al 2019 Phys. Rev. D 100 064003 (Preprint http://arxiv.org/abs/1901.07038)
  * [56] Ramos-Buades A, Husa S, Pratten G, Estellés H, García-Quirós C, Mateu-Lucena M, Colleoni M and Jaume R 2020 Phys. Rev. D 101 083015 (Preprint http://arxiv.org/abs/1909.11011)
  * [57] McWilliams S T 2019 Phys. Rev. Lett. 122(19) 191102 URL https://link.aps.org/doi/10.1103/PhysRevLett.122.191102
  * [58] Boetzel Y, Susobhanan A, Gopakumar A, Klein A and Jetzer P 2017 Phys. Rev. D 96 044011 (Preprint http://arxiv.org/abs/1707.02088)
  * [59] Fitzpatrick R 2012 An Introduction to Celestial Mechanics (Cambridge University Press)
  * [60] Colwell P 1993 Solving Kepler’s equation over three centuries
  * [61] Taylor S R, Huerta E A, Gair J R and McWilliams S T 2016 Astrophys. J. 817 70 (Preprint http://arxiv.org/abs/1505.06208)
  * [62] Klioner S A 2016 arXiv e-prints arXiv:1609.00915 (Preprint http://arxiv.org/abs/1609.00915)
  * [63] Arfken G B, Weber H J and Harris F E 2012 Mathematical methods for physicists: a comprehensive guide; 7th ed. (Burlington, MA: Elsevier) URL https://cds.cern.ch/record/1438219
  * [64] Einstein A, Infeld L and Hoffmann B 1938 Annals Math. 39 65–100
  * [65] Poisson E and Will C M 2014 Gravity
  * [66] Loutrel N, Liebersbach S, Yunes N and Cornish N 2019 Class. Quant. Grav. 36 01 (Preprint http://arxiv.org/abs/1801.09009)
  * [67] Damour T, Gopakumar A and Iyer B R 2004 Phys. Rev. D 70 064028 (Preprint http://arxiv.org/abs/gr-qc/0404128)
  * [68] Memmesheimer R M, Gopakumar A and Schaefer G 2004 Phys. Rev. D 70 104011 (Preprint http://arxiv.org/abs/gr-qc/0407049)
  * [69] Konigsdorffer C and Gopakumar A 2006 Phys. Rev. D 73 124012 (Preprint http://arxiv.org/abs/gr-qc/0603056)
  * [70] Klein A, Boetzel Y, Gopakumar A, Jetzer P and de Vittori L 2018 Phys. Rev. D 98 104043 (Preprint http://arxiv.org/abs/1801.08542)
  * [71] Moore B and Yunes N 2019 Class. Quant. Grav. 36 185003 (Preprint http://arxiv.org/abs/1903.05203)
  * [72] Boetzel Y, Mishra C K, Faye G, Gopakumar A and Iyer B R 2019 Phys. Rev. D 100 044018 (Preprint http://arxiv.org/abs/1904.11814)
  * [73] Hinder I, Herrmann F, Laguna P and Shoemaker D 2010 Phys. Rev. D 82 024033 (Preprint http://arxiv.org/abs/0806.1037)
  * [74] Goldberg J N, MacFarlane A J, Newman E T, Rohrlich F and Sudarshan E C G 1967 J. Math. Phys. 8 2155
  * [75] Buonanno A, Chen Y b and Vallisneri M 2003 Phys. Rev. D 67 104025 [Erratum: Phys.Rev.D 74, 029904 (2006)] (Preprint http://arxiv.org/abs/gr-qc/0211087)
  * [76] Huerta E A et al 2017 Phys. Rev. D 95 024038 (Preprint http://arxiv.org/abs/1609.05933)
  * [77] Berti E, Cardoso V, Gonzalez J A, Sperhake U, Hannam M, Husa S and Brügmann B 2007 Physical Review D 76 ISSN 1550-2368 URL http://dx.doi.org/10.1103/PhysRevD.76.064034
  * [78] Blanchet L, Qusailah M S S and Will C M 2005 The Astrophysical Journal 635 508?515 ISSN 1538-4357 URL http://dx.doi.org/10.1086/497332
  * [79] Bardeen J M, Press W H and Teukolsky S A 1972 Astrophys. J. 178 347
  * [80] Ori A and Thorne K S 2000 Physical Review D 62 ISSN 1089-4918 URL http://dx.doi.org/10.1103/PhysRevD.62.124022
  * [81] Buonanno A and Damour T 2000 Physical Review D 62 ISSN 1089-4918 URL http://dx.doi.org/10.1103/PhysRevD.62.064015
  * [82] Husa S and Winicour J 1999 Physical Review D 60 ISSN 1089-4918 URL http://dx.doi.org/10.1103/PhysRevD.60.084019
  * [83] Faraoni V and Prain A 2015 Understanding dynamical black hole apparent horizons (Preprint http://arxiv.org/abs/1511.07775)
  * [84] An X 2017 Emergence of Apparent Horizon in Gravitational Collapse (Preprint http://arxiv.org/abs/1703.00118)
  * [85] Price R H and Pullin J 1994 Phys. Rev. Lett. 72 3297–3300 (Preprint http://arxiv.org/abs/gr-qc/9402039)
  * [86] Pullin J 1999 Progress of Theoretical Physics Supplement 136 107?120 ISSN 0375-9687 URL http://dx.doi.org/10.1143/PTPS.136.107
  * [87] Le Tiec A and Blanchet L 2010 Class. Quant. Grav. 27 045008 (Preprint http://arxiv.org/abs/0910.4593)
  * [88] Khanna G and Price R H 2017 Phys. Rev. D 95 081501 (Preprint http://arxiv.org/abs/1609.00083)
  * [89] Berti E, Cardoso V and Will C M 2006 Physical Review D 73 ISSN 1550-2368 URL http://dx.doi.org/10.1103/PhysRevD.73.064030
  * [90] Sathyaprakash B S and Schutz B F 2009 Living Reviews in Relativity 12 2 (Preprint http://arxiv.org/abs/0903.0338)
  * [91] Nyadzani L and Razzaque S 2019 Coalescence Time of Compact Binary Systems (Preprint http://arxiv.org/abs/1905.06086)
  * [92] Zwick L, Capelo P R, Bortolas E, Mayer L and Amaro-Seoane P 2020 Mon. Not. Roy. Astron. Soc. 495 2321–2331 (Preprint http://arxiv.org/abs/1911.06024)
  * [93] Mishra C K, Kela A, Arun K G and Faye G 2016 Phys. Rev. D 93 084054 (Preprint http://arxiv.org/abs/1601.05588)
  * [94] Marchand T, Bernard L, Blanchet L and Faye G 2018 Phys. Rev. D 97 044023 (Preprint http://arxiv.org/abs/1707.09289)
  * [95] Blanchet L and Faye G 2019 Class. Quant. Grav. 36 085003 (Preprint http://arxiv.org/abs/1811.08966)
  * [96] Foffa S and Sturani R 2019 Phys. Rev. D 100 024047 (Preprint http://arxiv.org/abs/1903.05113)
  * [97] Henry Q, Faye G and Blanchet L 2020 Phys. Rev. D 102 044033 (Preprint http://arxiv.org/abs/2005.13367)
  * [98] Bini D, Damour T and Geralico A 2020 Phys. Rev. D 102 024062 (Preprint http://arxiv.org/abs/2003.11891)
  * [99] Bini D, Damour T and Geralico A 2020 Phys. Rev. D 102 024061 (Preprint http://arxiv.org/abs/2004.05407)
  * [100] Blümlein J, Maier A, Marquard P and Schäfer G 2021 Phys. Lett. B 816 136260 (Preprint http://arxiv.org/abs/2101.08630)
  * [101] Research W 2019 NDSolve https://reference.wolfram.com/language/ref/NDSolve.html accessed: 10-November-2021
  * [102] Inc W R Mathematica, Version 12.3.1 champaign, IL, 2021 URL https://www.wolfram.com/mathematica
  * [103] Arun K G, Blanchet L, Iyer B R and Sinha S 2009 Phys. Rev. D 80 124018 (Preprint http://arxiv.org/abs/0908.3854)
  * [104] Hwang H C, Ting Y S and Zakamska N L 2021 The eccentricity distribution of wide binaries and their individual measurements (Preprint http://arxiv.org/abs/2111.01789)
  * [105] Zevin M, Samsing J, Rodriguez C, Haster C J and Ramirez-Ruiz E 2019 Astrophys. J. 871 91 (Preprint http://arxiv.org/abs/1810.00901)
  * [106] Buskirk D and Hamilton M C B 2019 Eur. J. Phys. 40 025603 (Preprint http://arxiv.org/abs/1810.06160)
  * [107] Buskirk D P 2019 Marshall university marshall digital scholar URL https://mds.marshall.edu/cgi/viewcontent.cgi?article=2254&context=etd
  * [108] Christodoulou D 1970 Phys. Rev. Lett. 25(22) 1596–1597
  * [109] Echeverria F 1989 Gravitational-Wave Measurements of the Mass and Angular Momentum of a Black Hole General Relativity and Gravitation p 571


# Lackadaisical quantum walks on 2D grids with multiple marked vertices

Nikolajs Nahimovs and Raqueline A. M. Santos

(Center for Quantum Computer Science, Faculty of Computing, University of
Latvia
Raina bulv. 19, Riga, LV-1586, Latvia
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Lackadaisical quantum walk (LQW) is a quantum analog of a classical lazy walk,
where each vertex has a self-loop of weight $l$. For a regular
$\sqrt{N}\times\sqrt{N}$ 2D grid LQW can find a single marked vertex with
$O(1)$ probability in $O(\sqrt{N\log N})$ steps using $l=d/N$, where $d$ is
the degree of the vertices of the grid [12]. For multiple marked vertices,
however, $l=d/N$ is not optimal as the success probability decreases with the
increase of the number of marked vertices [13]. In this paper, we numerically
study search by LQW for different types of 2D grids – triangular, rectangular
and honeycomb – with multiple marked vertices. We show that in all cases the
weight $l=m\cdot d/N$, where $m$ is the number of marked vertices, still leads
to $O(1)$ success probability.

##  1 Introduction

Quantum walks are quantum counterparts of classical random walks [15].
Similarly to classical random walks, there are two types of quantum walks:
discrete-time quantum walks (DTQW), introduced by Aharonov et al. [3], and
continuous-time quantum walks (CTQW), introduced by Farhi et al. [10]. For the
discrete-time version, the step of the quantum walk is usually given by two
operators – coin and shift – which are applied repeatedly. The coin operator
acts on the internal state of the walker and rearranges the amplitudes of
going to adjacent vertices. The shift operator moves the walker between the
adjacent vertices.

Quantum walks have been useful for designing algorithms for a variety of
search problems [16]. To solve a search problem using quantum walks, we
introduce the notion of marked elements (vertices), corresponding to elements
of the search space that we want to find. We perform a quantum walk on the
search space with one transition rule at the unmarked vertices, and another
transition rule at the marked vertices. If this process is set up properly, it
leads to a quantum state in which the marked vertices have higher probability
than the unmarked ones. This method of search using quantum walks was first
introduced in [18] and has been used many times since then.

The concept of lackadaisical quantum walk (LQW), i.e. quantum walk with self-
loops, was first studied for DTQW on the one-dimensional line [14, 19]. Later
it was successfully applied to improve the DTQW based search on the complete
graph [21] and two-dimensional rectangular grid [22]. For the rectangular 2D
grid LQW gives $O(\sqrt{\log{N}})$ speed-up over the non-lackadaisical
algorithm of [5]111There are also other methods to achieve similar speed-up,
e.g. controlling the quantum walk using an ancilla qubit [20] or classically
searching the neighbourhood of the found vertex [4]..

The running time of the lackadaisical walk heavily depends on the weight of
the self-loop. For a regular 2D grid with a single marked vertex, the optimal
weight of the self-loop is $l=d/N$, where $d$ is the degree of the vertices of
the grid [12]. For multiple marked vertices, however, $l=d/N$ is not optimal
as the success probability decreases with the increase of the number of marked
vertices [13].

There are a few papers studying LQW search on the rectangular 2D grid with
multiple marked vertices. Saha et al. [17] showed that if $m$ marked vertices
are arranged in a $\sqrt{m}\times\sqrt{m}$ cluster, one should use
$l=\frac{4}{N\left(m+\lfloor{\frac{\sqrt{m}}{2}}\rfloor\right)}$. Giri and
Korepin [11] numerically found suitable choices of $l$ for which it is
possible to find up to 6 marked vertices in
$O\left(\sqrt{\frac{N}{m}\log\frac{N}{m}}\right)$ time steps with $O(1)$
success probability. Nahimovs [13] have demonstrated the existence of
exceptional configurations of marked vertices (i.e. configurations of marked
vertices for which the probability of finding a marked vertex does not grow
with the number of steps) and proposed two values of $l$ of the form
$\frac{4(m-O(m))}{N}$ which results in $O(1)$ success probability. Carvalho
_et al._ [9] numerically showed that for $N\gg m$ the success probability is
inversely proportional to the density of marked vertices and directly
proportional to the relative distance between the marked vertices.

In this paper, we numerically study search by LQW for different types of 2D
grids – triangular, rectangular and honeycomb – with multiple marked vertices.
We show that in all cases the weight $l=m\cdot d/N$, where $m$ is the number
of marked vertices, leads to $O(1)$ success probability. The results are
obtained from numerical simulations. The .NET code used to simulate the
quantum walk search algorithms is available on GitHub [6].

The paper is organized as follows. In Section 2, we define the lackadaisical
quantum walk and how we can do search on rectangular, triangular and honeycomb
2D grids. In Section 3, we find a suitable option for the self-loop weight and
analyse its behavior when searching an arbitrary set of marked vertices. And
we draw our conclusions in Section 4.

##  2 Lackadaisical Quantum walks on two-dimensional grids

Consider a two-dimensional $\\{$triangular, rectangular, honeycomb$\\}$ grid
of size $\sqrt{N}\times\sqrt{N}$ with periodic (torus-like) boundary
conditions. The vertices of the grid are labeled by the coordinates $(x,y)$
for $x,y\in\\{0,\dots,\sqrt{N}-1\\}$. The coordinates define a set of state
vectors, $\\{|x,y\rangle:x,y\in\\{0,\dots,\sqrt{N}-1\\}\\}$, which span the
$N$-dimensional Hilbert space ${\cal{H_{P}}}$ associated with the position.
Let $d$ be the degree of the grid. The lackadaisical quantum walk [23] has an
additional self-loop of weight $l$ in each vertex. Then, the Hilbert space
associated with the directions the walker can face is a $d+1$-dimensional
Hilbert space. We refer to it as the coin subspace ${\cal{H}_{C}}$. Therefore,
the Hilbert space of the lackadaisical quantum walk is
$\mathcal{H}={\cal{H_{P}}}\otimes{\cal{H_{C}}}=\mathbb{C}^{N}\otimes\mathbb{C}^{d+1}$.

The evolution of a state of the walk (without searching) is driven by the
unitary operator $U=S\cdot(I_{N}\otimes C)$, where $S$ is the flip-flop shift
operator [5] and $C$ is the coin operator, given by

$C=2|s_{c}\rangle\langle s_{c}|-I_{d+1},$ (1)

where

$|s_{c}\rangle=\frac{1}{\sqrt{d+l}}\sum_{i=0}^{d-1}|i\rangle+\sqrt{l}|\circlearrowleft\rangle.$

Notice that when $l=0$ we obtain the regular (non-lackadaisical) quantum walk.
Following we draw up the specifics for the lackadaisical quantum walk for each
type of grid.

###  2.1 The rectangular grid

The coin subspace of the walk (see Fig. 1) is a 5-dimensional Hilbert space
spanned by the set of states
$\\{|c\rangle:c\in\\{\uparrow,\downarrow,\leftarrow,\rightarrow,\circlearrowleft\\}\\}$.

Figure 1: The two-dimensional rectangular grid.

The shift operator acts as

$\displaystyle S|x,y,\uparrow\rangle$ $\displaystyle=$
$\displaystyle|x,y+1,\downarrow\rangle,$ $\displaystyle
S|x,y,\downarrow\rangle$ $\displaystyle=$
$\displaystyle|x,y-1,\uparrow\rangle,$ $\displaystyle S|x,y,\leftarrow\rangle$
$\displaystyle=$ $\displaystyle|x-1,y,\rightarrow\rangle,$ $\displaystyle
S|x,y,\rightarrow\rangle$ $\displaystyle=$
$\displaystyle|x+1,y,\leftarrow\rangle,$ $\displaystyle
S|x,y,\circlearrowleft\rangle$ $\displaystyle=$
$\displaystyle|x,y,\circlearrowleft\rangle.$

###  2.2 The triangular grid

The coin subspace of the walk is a 7-dimensional Hilbert space spanned by the
set of states
$\\{|c\rangle:c\in\\{\nwarrow,\nearrow,\leftarrow,\rightarrow,\swarrow,\searrow,\circlearrowleft\\}\\}$.

(a) (b)

(c)

Figure 2: (a) The two-dimensional triangular grid and (b)-(c) its mapping to
the rectangular grid.

There exists a simple mapping from the triangular to rectangular grid as shown
on Fig. 2. This allows us to label vertices of the grid by the coordinates
$(x,y)$ for $x,y\in\\{0,\dots,\sqrt{N}-1\\}$. The shift operator $S$ acts as

$\displaystyle S|x,y,\nwarrow\rangle$ $\displaystyle=$
$\displaystyle|x-1,y+1,\searrow\rangle,$ $\displaystyle S|x,y,\searrow\rangle$
$\displaystyle=$ $\displaystyle|x+1,y-1,\nwarrow\rangle,$ $\displaystyle
S|x,y,\leftarrow\rangle$ $\displaystyle=$
$\displaystyle|x-1,y,\rightarrow\rangle,$ $\displaystyle
S|x,y,\rightarrow\rangle$ $\displaystyle=$
$\displaystyle|x+1,y,\leftarrow\rangle,$ $\displaystyle S|x,y,\swarrow\rangle$
$\displaystyle=$ $\displaystyle|x,y-1,\nearrow\rangle,$ $\displaystyle
S|x,y,\nearrow\rangle$ $\displaystyle=$ $\displaystyle|x,y+1,\swarrow\rangle,$
$\displaystyle S|x,y,\circlearrowleft\rangle$ $\displaystyle=$
$\displaystyle|x,y,\circlearrowleft\rangle.$

###  2.3 The honeycomb grid

For the honeycomb grid there are two types of vertices, having either
$\\{\rightarrow,\nwarrow,\swarrow,\circlearrowleft\\}$ or
$\\{\leftarrow,\searrow,\nearrow,\circlearrowleft\\}$ directions. We can use a
$4$-dimensional coin space
$\\{\leftrightarrow,\mathrel{\text{\ooalign{$\swarrow$\cr$\nearrow$}}},\mathrel{\text{\ooalign{$\nwarrow$\cr$\searrow$}}},\circlearrowleft\\}$
which corresponds to either
$\\{\rightarrow,\nwarrow,\swarrow,\circlearrowleft\\}$ or
$\\{\leftarrow,\searrow,\nearrow,\circlearrowleft\\}$ depending on the type of
vertex.

(a) (b)

(c)

Figure 3: (a) The two-dimensional honeycomb grid and (b)-(c) its mapping to
the rectangular grid.

There exists a simple mapping from the honeycomb to the rectangular grid as
shown on Fig. 3. This allows us to label the vertices of the grid by the
coordinates $(x,y)$ for $x,y\in\\{0,\dots,\sqrt{N}-1\\}$. The shift operator
$S$ acts on the first type of vertices as

$\displaystyle S|x,y,\nwarrow\rangle$ $\displaystyle=$
$\displaystyle|x,y+1,\searrow\rangle,$ $\displaystyle S|x,y,\searrow\rangle$
$\displaystyle=$ $\displaystyle|x,y-1,\nwarrow\rangle,$ $\displaystyle
S|x,y,\leftarrow\rangle$ $\displaystyle=$
$\displaystyle|x-1,y,\rightarrow\rangle,$ $\displaystyle
S|x,y,\circlearrowleft\rangle$ $\displaystyle=$
$\displaystyle|x,y,\circlearrowleft\rangle$

and on the second type of vertices as

$\displaystyle S|x,y,\rightarrow\rangle$ $\displaystyle=$
$\displaystyle|x+1,y,\leftarrow\rangle,$ $\displaystyle S|x,y,\swarrow\rangle$
$\displaystyle=$ $\displaystyle|x,y-1,\nearrow\rangle,$ $\displaystyle
S|x,y,\nearrow\rangle$ $\displaystyle=$ $\displaystyle|x,y+1,\swarrow\rangle,$
$\displaystyle S|x,y,\circlearrowleft\rangle$ $\displaystyle=$
$\displaystyle|x,y,\circlearrowleft\rangle.$

###  2.4 Searching

In order to search for marked vertices we extend the step of the algorithm,
making it

$U^{\prime}=U\cdot(Q\otimes I_{d+1}),$

where $Q$ is the query transformation which flips the sign at marked vertices,
irrespective of the coin state. The system starts in

$|\psi(0)\rangle=\frac{1}{\sqrt{N}}\sum_{x,y=0}^{\sqrt{N}-1}|x,y\rangle\otimes|s_{c}\rangle,$
(2)

which is the uniform superposition over vertices (but not directions). Note
that $|\psi(0)\rangle$ is a 1-eigenvector of $U$ but not of $U^{\prime}$. The
state of the system after $t$ steps is $|\psi(t)\rangle=U^{\prime
t}|\psi(0)\rangle$. The probability of finding a marked vertex at time $t$ is
given by

$p_{M}(t)=\langle\psi(t)|\left(\sum_{(x,y)\in M}|x,y\rangle\langle x,y|\otimes
I_{d+1}\right)|\psi(t)\rangle,$ (3)

where $M$ is the set of marked vertices. If there are marked vertices, the
state of the algorithm starts to deviate from $|\psi(0)\rangle$. For the
regular quantum walk ($l=0$), in the case of a single marked vertex, after
$O(\sqrt{N\log{N}})$ steps the inner product $\langle\psi(t)|\psi(0)\rangle$
becomes close to $0$. If the state is measured at this moment, the probability
of finding a marked vertex is $O(1/\log{N})$ [5, 2, 1]. By using amplitude
amplification [8] we obtain the total running time of $O(\sqrt{N}\log{N})$
steps. For the lackadaisical quantum walk, we can find a single marked vertex
with $O(1)$ probability in $O(\sqrt{N\log N})$ steps using $l=d/N$ [12].

##  3 Finding a suitable weight for multiple marked vertices

The running time of the search by LQW depends on the weight of the self-loop
$l$. For a regular 2D grid with a single marked vertex, the optimal weight of
the self-loop is $l=d/N$ [12]. For multiple marked vertices, $l=d/N$ is not
optimal.

In this section we analyse the effect of different values of $l$ for the three
types of 2D grids: triangular, rectangular and honeycomb. To start, let us
plot the probability of finding a marked vertex for different values of the
self-loop weight $l$.

(a) (b) (c)

Figure 4: The probability of finding a marked vertex for different values of
$l$ for the (a) triangular, (b) rectangular and (c) honeycomb grids of size
$24\times 24$ with set of marked vertices
$\\{(0,0),(0,2),(0,4),(0,6),(0,8)\\}$. The solid blue curve is $l=0$; the
dotted yellow curve is $l=d/N$; the dashed green curve is the optimal value
for $l$; the dot-dashed red curve is $l=d\cdot m/N$; and the dot-dot-dashed
purple curve is $l=0.1$.

Figure 4 shows the triangular (4a), rectangular (4b) and honeycomb (4c) grids
of size $N=24\times 24$ with $m=5$ marked vertices
$\\{(0,0),(0,2),(0,4),(0,6),(0,8)\\}$. We can observe the same behavior for
all grids. At $t=0$, the success probability is equal to $m/N=5/576\approx
0.0087$. The solid blue curve $l=0$ depicts the evolution of the regular (non-
lackadaisical) quantum walk. If we increase the value of $l$ the success
probability starts to grow. The dotted yellow curve depicts $l=d/N$ – the
optimal value for a single marked vertex. Since we have more than one marked
vertex this value of $l$ is not optimal. The green dashed curve shows the
success probability, corresponding to $l=l_{opt}\approx 0.0490,0.0339,0.0207$
for the triangular, rectangular and honeycomb grids, respectively. The dot-
dashed red curve $l=d\cdot m/N$ is very close to the optimal curve, specially
for the rectangular grid. By further increasing the value of $l$, the success
probability starts to drop, as depicted by the purple dot-dot-dashed curve
$l=0.1$.

Figure 5 compares three values for the self-loop weight $l=a\cdot d/N$,
$a\in\\{m,m-\sqrt{m},1\\}$. Fig. 5a, 5b and 5c plot the maximum success
probability (first peak) for the triangular, rectangular and honeycomb grids,
respectively. The grids size is $N=100\times 100$. The set of marked vertices
is $M_{m}=\\{(0,2i)|i\in[0,m-1]\\}$ for $m\in[1,48]$. The dotted green curves
for $a=1$ indicates a rapid decrease of the success probability as $m$
increases. It confirms that $l=d/N$ is not a good choice for improving the
search algorithm. The other two values, $a=m$ and $a=m-\sqrt{m}$, look more
promising. They were first proposed in [13] for the rectangular grid. We can
observe that the solid blue curves for $a=m$ and dashed yellow curves for
$a=m-\sqrt{m}$ get close to each other as $m$ increases. Comparing both
values, the success probability for $a=m-\sqrt{m}$ is considerably smaller for
small values of $m$ and it gets slightly bigger for big values of $m$. Overall
$a=m$ behaves better and it can be a suitable choice for the self-loop weight
in order to improve the search algorithm. In the next section, we will analyse
the behavior of the self-loop weight $l=m\cdot d/N$ when searching for an
arbitrary set of marked vertices.

(a) (b) (c)

Figure 5: Comparison of the self-loop weights $l=a\cdot d/N$, where
$a\in\\{m,m-\sqrt{m},1\\}$, for $100\times 100$ grid with set of marked
vertices $M_{m}=\\{(0,2i)|i\in[0,m-1]\\}$ for $m\in[1,48]$. The maximum
success probability is shown for the (a) triangular, (b) rectangular and (c)
honeycomb grids. The solid blue curves depicts $a=m$; the dashed yellow
curves, $a=m-\sqrt{m}$; and the dotted green curves, $a=1$.

###  3.1 Analysis for the self-loop weight $l=m\cdot d/N$

The following plots show the average over 100 executions of the search
algorithm, considering randomly generated sets of marked vertices for
$m\in[1,1000]$ for a grid of size $100\times 100$. Figure 6 depicts the
average maximum success probability as a function of the number of marked
vertices $m$. For all three types of grid the average success probability is
higher and very close to 1 for $m=1$ and it is inversely proportional to the
number of marked vertices (as noticed before in [9] for the rectangular grid).
For the shown range of marked vertices, the success probability stays above
0.45. The standard deviation is very small, which is shown by the confidence
interval depicted by the light blue area.

(a) (b) (c)

Figure 6: The average maximum success probability over 100 runs for the (a)
triangular, (b) rectangular and (c) honeycomb grids. The self-loop weight is
$l=m\cdot d/N$. The grids size is $N=100\times 100$. The light blue area
depicts the confidence band. The dotted light green and red curves are the
maximum and minimum values, respectively, over the 100 runs. The sets of
marked vertices with sizes ranging from 1 to 1000 are randomly generated.

It is worth to mention that in our experiments sets of marked vertices were
generated at random and we didn’t prevent them to form an exceptional
configuration222It is known, that lackadaisical quantum walks have exceptional
configurations [13] of marked vertices. In such cases, the probability of
finding a marked vertex stays close to the initial probability during the
evolution.. This, however, did not affect the results. This can be seen from
the values of the minimum probability (dotted light red curves) which are far
from the initial probability. Note, that it is still possible that some
smaller subsets of marked vertices have formed exceptional configurations.

The average number of steps (the time to reach the first peak of the success
probability) required by the algorithm is shown by Figs. 7a, 7b and 7c for the
rectangular, triangular and honeycomb grids, respectively. The plots are in
log-log scale. The dashed yellow curve represents the function
$\sqrt{\frac{N}{m}\log{\frac{N}{m}}}$. We can see that the curves for the
number of steps and the dashed yellow curve have almost the same inclination.
For the triangular and rectangular grids, the dashed yellow curve is almost a
perfect fit. The same behavior can be seen when we scale the size of the
grids. This means that the number of steps required by the search algorithm is
$O\left(\sqrt{\frac{N}{m}\log{\frac{N}{m}}}\right)$.

(a) (b) (c)

Figure 7: The average number of steps (solid blue curves) over 100 cases for
the (a) triangular, (b) rectangular and (c) honeycomb grids in log-log scale.
The self-loop weight is $l=m\cdot d/N$. The grids size is $N=100\times 100$.
The light blue area depicts the confidence band. The dashed yellow curve is
$\sqrt{\frac{N}{m}\log\frac{N}{m}}$. The sets of marked vertices with sizes
ranging from 1 to 1000 are randomly generated.

Let us see what happens when we increase the number of marked vertices even
further. Consider $20\%$ of vertices marked, that is, $m=N/5$. Figure 8 shows
the average maximum success probability over 100 cases for the 3 types of grid
with different sizes. We can see that the values stays above 0.4 for all
grids. The honeycomb grid is the one with lower values, followed by the
rectangular and then the triangular grid. The standard deviation is very small
and it can be barely seen by the error bars.

Figure 8: The average maximum success probability over 100 cases for different
sizes of the triangular, rectangular and honeycomb grids. The self-loop weight
is $l=m\cdot d/N$. The sets of marked vertices with size $m=N/5$ are randomly
generated.

In conclusion, since we know that the success probability is above some
threshold, we can search with $O(1)$ probability in
$O\left(\sqrt{\frac{N}{m}\log{\frac{N}{m}}}\right)$ time.

##  4 Conclusions

We have studied the search for multiple marked vertices by lackadaisical
quantum walks on triangular, rectangular and honeycomb 2D grids. Previous to
our work, a few researchers [17, 11, 13, 9] have studied the case of having
multiple marked vertices. However, the previous research was focused on the
rectangular grid and covered only a small range and specific sets of marked
vertices. In this paper, we presented numerical results for different types of
grids and a broad range and arbitrary sets of marked vertices. We have
demonstrated that by using lackadaisical quantum walks with self-loop weight
of $l=m\cdot d/N$, we can find a marked vertices on triangular, rectangular
and honeycomb 2D grids with $O(1)$ probability in
$O\left(\sqrt{\frac{N}{m}\log\frac{N}{m}}\right)$ steps. The presented results
are based on numerical simulations and it would be nice to obtain analytical
results using the new method developed in [7].

## Acknowledgements

This work was supported by ERDF project number 1.1.1.2/VIAA/1/16/002 and by
the QuantERA ERA-NET Cofund in Quantum Technologies implemented within the
European Union’s Horizon 2020 Programme (QuantAlgo project).

## References

  * [1] G. Abal, R. Donangelo, M. Forets, and R. Portugal.  Spatial quantum search in a triangular network.  Mathematical Structures in Computer Science, 22:521–531, Jun 2012\.
  * [2] G. Abal, R. Donangelo, F. L. Marquezino, and R. Portugal.  Spatial search on a honeycomb network.  Mathematical Structures in Computer Science, 15:999–1009, Dec 2010\.
  * [3] Y. Aharonov, L. Davidovich, and N. Zagury.  Quantum random walks.  Physical Review A, 48(2):1687–1690, 1993.
  * [4] A. Ambainis, A. Bačkurs, N. Nahimovs, R. Ozols, and A. Rivosh.  Search by quantum walks on two-dimensional grid without amplitude amplification.  In Proceedings of the 7th Annual Conference on Theory of Quantum Computation, Communication, and Cryptography (TQC2012), page 87–97. Springer, Tokyo, Japan, 2013.
  * [5] A. Ambainis, J. Kempe, and A. Rivosh.  Coins make quantum walks faster.  In Proceedings of the 16th ACM-SIAM Symposium on Discrete Algorithms, pages 1099–1108, 2005.
  * [6] CQCS at University of Latvia.  Quantum walk simulator.  URL: https://github.com/CQCS-LU/Quantum-walk-simulator-NET.
  * [7] G. A. Bezerra, P. H. G. Lugão, and R. Portugal.  Quantum walk-based search algorithms with multiple marked vertices, 2021\.  arXiv:2103.12878.
  * [8] Gilles Brassard, Peter Høyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation, pages 53–74.  2002\.
  * [9] Jonathan H. A. de Carvalho, Luciano S. de Souza, Fernando M. de Paula Neto, and Tiago A. E. Ferreira.  Impacts of Multiple Solutions on the Lackadaisical Quantum Walk Search Algorithm, pages 122–135.  Springer International Publishing, 2020.
  * [10] E. Farhi and S. Gutmann.  Quantum computation and decision trees.  Physical Review A, 58:915–928, 1998.
  * [11] Pulak Ranjan Giri and Vladimir Korepin.  Lackadaisical quantum walk for spatial search.  Modern Physics Letters A, 35(08):2050043, 2020.
  * [12] Peter Høyer and Zhan Yu.  Analysis of lackadaisical quantum walks.  Quantum Inf. Comput., 20(13&14):1137–1152, 2020.
  * [13] N. Nahimovs.  Lackadaisical quantum walks with multiple marked vertices.  In Proceedings of SOFSEM 2019, volume 11376, pages 368–378, 2019\.
  * [14] Inui Norio, Konno Norio, and Segawa Etsuo.  One-dimensional three-state quantum walk.  Physical Review E. Statistical Nonlinear and Soft Matter Physics, 72(5 Pt 2):168–191, 2005.
  * [15] R. Portugal.  Quantum walks and search algorithms.  Springer, New York, 2013.
  * [16] D. Reitzner, D. Nagaj, and V. Buzek.  Quantum walks.  In Acta Physica Slovaca 61, volume 6, pages 603–725, 2011.  arxiv.org/abs/1207.7283.
  * [17] A. Saha, R. Majumdar, D. Saha, A. Chakrabarti, and S. Sur-Kolay.  Search of clustered marked states with lackadaisical quantum walks.  ArXiv e-prints, 2018.  arXiv:1804.01446.
  * [18] N. Shenvi, J. Kempe, and K. B. Whaley.  A quantum random walk search algorithm.  Physical Review A, 67(052307), 2003.
  * [19] M. Stefanak, I. Bezdekova, and I Jex.  Limit distributions of three-state quantum walks: the role of coin eigenstates.  Physical Review A, 90(1):124–129, 2014.
  * [20] A. Tulsi.  Faster quantum-walk algorithm for the two-dimensional spatial search.  Physical Review A, 78(012310), 2008.
  * [21] T. G. Wong.  Grover search with lackadaisical quantum walks.  Journal of Physics A Mathematical General, 48, 2015.
  * [22] T. G. Wong.  Faster search by lackadaisical quantum walk.  Quantum Information Processing, 17, 2018.
  * [23] T. G. Wong and A. Ambainis.  Quantum search with multiple walk steps per oracle query.  Physical Review A, 92:10, Aug 2015.