# A Sea of Black Holes: Characterizing the LISA Signature for Stellar-Origin
Black Hole Binaries

Krystal Ruiz-Rocha  Department of Physics and Astronomy, Vanderbilt
University, Nashville, TN 37235, USA  Kelly Holley-Bockelmann  Department of
Physics and Astronomy, Vanderbilt University, Nashville, TN 37235, USA
Department of Life and Physical Sciences, Fisk University  Karan Jani
Department of Physics and Astronomy, Vanderbilt University, Nashville, TN
37235, USA  Michela Mapelli  Institut für Theoretische Astrophysik, ZAH,
Universität Heidelberg, Albert-Ueberle-Straße 2, D-69120,Heidelberg, Germany
Physics and Astronomy Department Galileo Galilei, University of Padova, Vicolo
dell’Osservatorio 3, I–35122, Padova, Italy  INFN - Padova, Via Marzolo 8,
I–35131 Padova, Italy  Samuel Dunham  TAPIR, Mailcode 350-17, California
Institute of Technology, 1200 E California Blvd, Pasadena, CA 91125, USA
William Gabella  Department of Physics and Astronomy, Vanderbilt University,
Nashville, TN 37235, USA

###### Abstract

Observations by the LIGO, Virgo and KAGRA (LVK) detectors have provided new
insights in the demographics of stellar-origin black hole binaries (sBHB). A
few years before gravitational-wave signals from sBHB mergers are recorded in
the LVK detectors, their early coalescence will leave a unique signature in
the ESA/NASA mission Laser Interferometer Space Antenna (LISA). Multiband
observations of sBHB sources between LISA and LVK detectors opens an
unprecedented opportunity to investigate the astrophysical environment and
multi-messenger early-alerts. In this study, we report the sBHB sources that
will be present in the LISA data derived directly from the hydrodynamic
cosmological simulation Illustris. By surveying snapshots across cosmological
volume, metallicity and look-back time, we find that about tens to thousand
sBHB candidates will be present in the LISA data for various combinations of
mission lifetime. For estimates consistent with the LVK rates, we find that
only 20 sBHBs across Illustris snapshots will be detected with significant
confidence for a 10-year LISA mission, while a 4-year LISA mission would
detect only 2 sBHBs. Our work paves the way for creating LISA mock data and
bench marking LISA detection pipelines directly using cosmological
simulations.

gravitational waves—stars:black holes—multimessenger astronomy—LISA

##  1 Introduction

The Laser Interferometer Gravitational-Wave Observatory (LIGO) made the first
direct detection of gravitational waves in 2015, revealing a merger of black
holes with masses 36 and 29 $M_{\odot}$ (Abbott et al., 2016). These black
holes, dubbed stellar origin binary black holes (sBHB) by the gravitational
wave community, are much more massive than the roughly 10 $M_{\odot}$ black
holes in accreting binaries that had been detected electromagnetically (van
den Heuvel, 2019). This discovery expanded our understanding of black hole
populations, providing insight into their formation and evolution (Abbott et
al., 2016a, b). However, the formation channel of these binaries is still
under debate. The current favored mechanisms involve formation via: a)
isolated binaries in the galactic field (Bethe & Brown, 1998; Belczynski et
al., 2002; Banerjee et al., 2010; Mennekens & Vanbeveren, 2014; Belczynski et
al., 2016; de Mink & Mandel, 2016; Marchant et al., 2016; Giacobbo et al.,
2018; Kruckow et al., 2018) ; b) dynamical interactions within dense stellar
environments (Portegies Zwart & McMillan, 2000; Banerjee et al., 2010; Ziosi
et al., 2014; Rodriguez et al., 2015, 2016; Mapelli, 2016; Banerjee, 2017;
Askar et al., 2017; Rodriguez et al., 2021); and c) binaries within gaseous
AGN accretion disks (McKernan et al., 2018; Stone et al., 2017; Bartos et al.,
2017; Ford & McKernan, 2022).

sBHBs are so massive that they are a new source class for the Laser
Interferometer Space Antenna (LISA), an ESA/NASA space-based gravitational
wave observatory set to launch in 2035 (Amaro-Seoane et al., 2017a; Colpi et
al., 2024). LISA is sensitive from around $10^{-4}$ to 0.1 Hz, a frequency
range that captures the inspiral phase of the sBHB population. In the LISA
band, low frequency sBHBs are relatively widely-separated, take millions of
years to merge, and may generate an unresolved stochastic background (Babak et
al., 2023). The higher frequency sources are more tightly-bound, sweeping
through the LISA band in the last phase of inspiral months to days before
merging within ground-based gravitational wave detectors, making them a prime
multiband source (Sesana, 2016; Cutler et al., 2019; Jani et al., 2020).

Multiband observations alone enable a more accurate measurement of masses,
spins, eccentricity, and distance (Vitale, 2016; Randall & Xianyu, 2021;
Toubiana et al., 2021; Ranjan et al., 2024), but their unique power is in
providing better spatial resolution and enough lead time before merger to
conduct rigorous observational campaigns to identify electromagnetic
counterparts (Baker et al., 2019; Lamberts et al., 2019; Korol et al., 2017;
Digman & Cornish, 2023). Multiband data is vital to differentiate between sBHB
formation channels as well as to constrain their final evolutionary sequence
(Breivik et al., 2016; Baker et al., 2019; Amaro-Seoane et al., 2017b).

Figure 1: Characteristic strain versus frequency for sBHBs formed within the
oldest snapshot of the Illustris-1 simulation at redshift z=13.4 (Mapelli et
al., 2017). We only show a subset of these sBHBs that eventually merged under
$z\leq 0.5$. The strain from each sBHB is colored by the total binary mass in
the source frame. The left curve line is LISA’s sensitivity curve, the right
curve is LIGO’s sensitivity curve.

To investigate the possibility that stellar origin binaries arise in an
isolated binary stellar population, Mapelli et al. (2017) used the
cosmological hydrodynamic simulation suite Illustris (Vogelsberger et al.,
2014a) to predict the astrophysical and LIGO-observable rate of sBHB mergers,
given a range of binary stellar population models. Illustris provided a self-
consistent time-evolving map of star-formation rate, mass, and metallicity to
be used as input to a binary population synthesis code, making it possible to
paint a stellar population onto the simulation, let it evolve, and track the
formation and ultimate merger of sBHBs.

Our work extends this study, identifying those sBHB mergers within the
Illustris simulation volume that are observable with LISA. Formally, the time
to merge for sBHBs in the LISA band can be from months to millions of years,
but we narrowed our investigation to those sBHBs that transition from LISA to
ground-based detectors within the lifetime of the LISA mission; these sources
are the loudest in the LISA band and are therefore the most promising for
multiband and multi-messenger observations (Gerosa et al., 2019). Figure 1 is
a subset of sBHB mergers within the Mapelli et al. (2017) database that may be
seen by LISA and merge within LIGO. These binaries are chosen from the
snapshot at the highest redshift(z = 13.4), because it contains sBHBs that
merge within a wide variety of lookback times, showing the robustness of
Mapelli et al. (2017)’s database for studying multi-messenger astronomy.

The paper proceeds as follows: Section 2 describes our methods, including a
brief review of the simulation, the stellar population input from Mapelli et
al. (2017). In addition, section 2.1 discusses the signal-to-noise
calculations for our binary database. In section 3 we discuss the results,
including the multiband rates (section 3.1), and number of detectable
multiband sources (section 3.2). Finally, the conclusion and discussion are
within section 4 where we summarize our results, point out caveats, and
identify areas for future work.

##  2 Methods

Illustris is a set of cosmological hydrodynamic N-body simulations run using
the AREPO moving mesh code (Vogelsberger et al., 2014a; Springel, 2010). To
model galaxy formation and evolution in a cosmological context, the simulation
self-consistently treats dark matter, supermassive black holes, gas, and
stars, including prescriptions for gas cooling, star formation, metal
enrichment, supernovae and associated feedback, as well as massive black hole
formation and evolution (including gas accretion, mergers, and several modes
of feedback). Illustris is a cosmological volume of $106.5^{3}$ Mpc3
initialized in a Wilkinson Microwave Anisotropy Probe (WMAP-9) cosmology with
the following parameters (Hinshaw et al., 2013): $\Omega_{M}=0.2726$,
$\Omega_{\Lambda}=0.7274$, $\Omega_{b}=0.0456$, h = $0.704$, and $H_{0}=70.4$
kms-1Mpc-1, where $\Omega$ represents the mean density relative to the
critical density of the Universe, and $M$, ${\Lambda}$, and $b$ refer to
matter, the cosmological constant, and baryons, respectively. $h$ is the
dimensionless Hubble parameter (Croton, 2013), and $H_{0}$ the Hubble
constant. The Illustris suite is composed of 6 simulations: Illustris-(1,2,3),
which are the hydrodynamical simulations modeling galaxy formation and include
all the star formation, gas and massive black hole physics described above,
and Illustris-DM-(1,2,3), which are the dark matter-only variant. We use
Illustris-1 for its higher mass resolution; the Illustris-1 simulation
contains $1.6\times 10^{12}$ particles over all snapshots, with a dark matter
mass resolution of $6.3\times 10^{6}M_{\odot}$, a baryonic matter mass
resolution of $1.3\times 10^{6}M_{\odot}$ and a softening length of $\sim 710$
pc for stars and supermassive black holes (SMBHs) (Vogelsberger et al., 2014a,
b).

Although the Illustris simulation suite is state-of-the-art, it cannot resolve
individual stars and stellar remnants. Therefore, Mapelli et al. (2017) maps a
population of sBHB, each stellar particle in snapshots spanning redshifts 0–16
using an updated version of the public binary stellar population synthesis
code, Binary Star Evolution (BSE) (Hurley et al., 2000, 2002) called MOBSE
(Giacobbo et al., 2018). For more detail about the procedure to seed Illustris
with sBHB, please see Mapelli et al. (2017); we briefly describe the process
here. Mapelli et al. (2017) constructed a library of 72 distinct binary
stellar population models, each ’book’ in the library consisting of millions
of sBHBs representing a particular metallicity and a set of assumptions about
the less well-constrained binary-evolution physics, such as mass loss in the
common-envelope phase.

Figure 2: Rate versus redshift for the fiducial model D and DK. The dotted
grey line and contours are a reproduction of LIGO’S O3 rate from Abbott et al.
(2023a)

Using this library, a sBHB distribution is assigned to Illustris-1 star
particles. The metallicity of a star particle determines which stellar
population ’book’ to use, and the mass of the star particle determines how
many sBHB to select. Each sBHB is characterized by its masses, separation, and
eccentricity, which is used to calculate the binary merger timescale. By
following only the sBHB that merge within a Hubble time, this earlier work
resulted in a simulated census of the astrophysical merger rates for these
sources, as well as estimates of the number of binaries and mass distribution
of sBHB mergers that LIGO will see for various stellar population models, as
shown in Figure 6 of Mapelli et al. (2017).

Note that for our paper, we adopt models D and DK as our fiducials; these
assume a delayed supernova model, a treatment of Hertzsprung gap (HG) donors,
and natal kicks for compact binaries (Mapelli et al., 2017). The difference
between model D and DK is that D employs a method described in Fryer et al.
(2012) to rescale natal kicks for black holes to take into account the
fraction of the stellar envelope that can fall back onto the compact object
during the explosion, $f_{b}$, which can range from 0 to 1 (Fryer et al.,
2012; Spera et al., 2015), as shown by equation 1

$v_{BH}=v_{NS}(1-f_{b}),$ (1)

where $v_{NS}$ is the velocity distribution in (Hobbs et al., 2005). Model DK
applies the velocity distribution found by Hobbs et al. (2005) for natal kicks
of neutron stars to the black holes.

Figure 3: Heatmap of BHB mass distribution for first BH mass($m_{1}$) and
second BH mass($m_{2}$) for the Mapelli post-processing files restricted to a
sBHBs merging below redshift 0.5. The original database did not have $m_{1}$
always greater than $m_{2}$. However within this figure we set $m_{1}$ to be
greater than $m_{2}$ for all sBHBs. The colors represent the log of the counts
of the sBHBs in each bin.

The sBHB merger rate in model D is $181$ $\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$
at a redshift of 0.2. At the time of publication, this rate was the most
consistent with O1 LIGO observations. However, LIGO data from its observing
run 3 (O3) places the sBHB merger rate between 17.9 and 44
$\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$ at a redshift of 0.2 (Abbott et al.,
2023a). This filters into our results as an overestimate of the LISA and LIGO
multiband rate when we use model D. In contrast, model DK offers a much closer
match to O3 rates with a sBHB merger rate of $29$
$\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$ at redshift 0.2. This similarity is shown
in Figure 2 which plots the merger rates for model D, DK, and LIGO’S rate for
O3 against redshift. The dashed line represents model D’s rate and the lower
line is model DK’s rate. LIGO’s rate is shown by the dotted gray line with the
inner contours containing 50% of the data and the outer contours containing
90% of the data. Model DK’s rate consistently falls within the 50% contours
across the full redshift range and aligns more closely to LIGO’s rate for
$z>0.3$. We note that the star formation rate history in the Illustris-1
simulation is overall in good agreement with observations (Madau & Fragos,
2017), except for an overestimate of $\sim{40\%}$ at low redshift
($z\lesssim{0.2}$). However, this overestimate has a negligible impact on the
merger rate density of sBHBs (Mapelli et al., 2017), because only stars with
nearly solar metallicity form at such a low redshift: the efficiency of BHB
mergers drops by about three orders of magnitude if the metallicity $Z$ is
$>0.2$ Z⊙ (Giacobbo et al., 2018).

###  2.1 Calculating the Signal-to-Noise Ratio in the LISA Band

We began with 108 snapshots ranging from redshift $0.0072$ to $13.6$, each
consisting of a database of the sBHB formed at that redshift. This database
contains information on: 1. the progenitor stellar particle ID from Illustis;
2. metallicity; 3. black hole component masses; 4. formation redshift; 5. sBHB
merger time.

Since each snapshot contains only those sBHB that formed within a particular
snapshot time interval, the varying sBHB merger time implies that a sBHB at a
previous snapshot may still exist at a given redshift, but will be missing in
the snapshot; indeed, the sBHBs seeded in a given snapshot will merge in the
future over a wide range of timescales. It is worth reiterating that the
database only contains sBHBs that merge within a Hubble time. Consequently,
sBHBs in the simulated population that do not merge by the present day are
missing from the current database. We leave the analysis of the quasi-static
sBHB population for a future paper.

Figure 4: Heatmap of the total mass and luminosity distance colored by the
average SNR in each grid for a LISA mission lifetime of 10 years. Figure 5:
Multiband rates calculated for our sBHBs database using equation 3 for both
Model D and Model DK for a 4-year LISA mission, shown with purple contours,
and a 10-year LISA mission shown with the teal contours.

To identify sBHB mergers that will be observable by both LISA and LIGO, we
first restrict the database to those sBHBs that merge within a horizon
redshift of 0.5; this decision greatly reduced the number of sBHBs in the
Illustris volume from $7.7\times 10^{9}$ to $1.19\times 10^{9}$. Our choice is
conservative and encompasses the expected horizon distance for sBHBs (see
Figure 2 of Jani et al. (2020)). Figure 3 displays the black hole mass
distribution for those systems that merge within a Hubble time and merge below
a redshift of 0.5 in model D. We also forced the $m_{1}$ to always be greater
than $m_{2}$, which was not the case in our database. In these models, the
most common sBHB merger arises between black holes with commensurate masses of
roughly 15–25 $M_{\odot}$. We caution that neither model D nor DK contain
sBHBs with component masses $m_{\rm BH}>43M_{\odot}$. LIGO has certainly
observed black holes more massive than the upper limit of these models (Abbott
et al., 2019, 2021, 2023b, 2024). It may well be that the high mass end of the
observed BHB merger mass spectrum represents additional processes, such as
hierarchical merging, that may be at play and are therefore not well-
represented by BSE (Hurley et al., 2000, 2002; Giacobbo et al., 2018). We
leave a multiband treatment of the higher mass black holes observed by LIGO to
a future paper.

Using the redshift-limited database, we extract the masses and lookback time
for each black hole binary merger and determine the gravitational wave
frequency, $f_{\rm obs}$, of these binaries at a time $T_{\rm obs}$ before the
merger, using:

$f_{\rm obs}=\frac{c^{3}}{8\pi
GM_{c}}\left(\frac{5GM_{c}}{c^{3}}\right)^{3/8}\frac{1}{T_{\rm obs}^{3/8}},$
(2)

where c is the speed of light, G is the gravitational constant, $M_{c}$ is the
chirp mass, and $T_{\rm obs}$ is the time until the BHB merges in seconds
(Chan et al., 2018; Bassan, 2014). We adopt 4 and 10 years for $T_{\rm obs}$,
the projected duration of science operations for the LISA mission. As outlined
in Gerosa et al. (2019), this choice produces a louder SNR in the LISA band
while enabling a relatively contemporaneous multiband source. (Gerosa et al.
(2019) considered a $T_{\rm wait}$, the time between detecting an event with a
space-based and a ground-based detector. In our work, $T_{\rm obs}$ is the
same as $T_{\rm wait}$.) The masses in equation 2 are in the detector frame,
and are converted from the source frame found in the database to detector
frame masses using: $m_{\rm detector}=m_{\rm source}(1+z)$, where z is the
sBHB merger redshift.

##  3 Results

###  3.1 Multiband Rate for detectable sBHBs

With this information we can then calculate the signal-to-noise ratio (SNR).
For this, we use the python LISA Sensitivity package given an observation
duration and cosmology (Robson et al., 2019). This code models LISA’s
sensitivity to laser shot noise, acceleration noise, and test-mass force noise
over a mission duration, and calculates the characteristic strain, and SNR of
a binary given the component masses, distance, and an optional sky location
and inclination. Note that unless specified, the SNR is averaged over
inclination angle, polarization, and sky position. Additionally, since the
sBHBs did not include spin, the gravitational wave signal is approximated by a
phenomenological waveform model for non-spinning binaries, referred to as
PhenomA (Ajith et al., 2007). In the database, the LISA SNR for a observation
run of 10 years technically ranged from 0 to $1\times 10^{7}$, with the
anomalously high SNR populated by a sBHB at a luminosity distance of less than
a kpc. Figure 4 maps the average SNR as a function of total mass and
luminosity distance assuming a mission duration of 10 years. In this figure we
see, as expected, that nearby high mass binaries yield the highest SNR. For
example sBHB system with $M_{\rm tot}=80\,M_{\odot}$ will visible up to 200
Mpc in LISA with SNR $\gtrsim 20$. Such high-mass binaries will be visible in
LISA out to ${\sim}1$ Gpc with an average SNR ${\sim}6$.

The luminosity distance for high-mass sBHB in Figure 4 at a SNR threshold of 8
is consistent with the multiband detection radius reported in Jani et al.
(2020).

Figure 6: The expected number of multiband sBHBs calculated for different SNR
thresholds. The left panel assumes model D and a 4 or 10-year LISA mission,
while the right panel model assumes model DK, which better approximates the O3
LIGO event rate.

Now that we have the SNR for each sBHB in hand, we can determine the multiband
event rates and the expected number of detections for LIGO and LISA as a
function of SNR.

We calculate the expected astrophysical rate of a population above a certain
SNR threshold using equation 3 from (Chen et al., 2022):

$\frac{dN}{dzdt}=\frac{d^{2}n(z)}{dzdV_{c}}\frac{dz}{dt}\frac{dV_{c}}{dz}\frac{1}{1+z},$
(3)

where $V_{c}$ is the co-moving volume, $n(z)$ is the number of sBHBs at a
redshift z, $\frac{1}{1+z}$ converts $\frac{dz}{dt}$ to the observer frame,
and $\frac{d^{2}n(z)}{dzdV_{c}}$ is the rate found in the simulation volume
above a certain SNR threshold. This can be found directly from the simulation,
expressed by:

$\frac{d^{2}n(z)}{dzdV_{c}}=\frac{N(z)}{\Delta zV_{sim}},$ (4)

where $N(z)$ is the number of sBHBs merging within a redshift bin, $\Delta z$
is the width of the redshift bin and $V_{\rm sim}$ is the simulation volume,
which in our case is $(106.5{\rm Mpc})^{3}$.

The multiband rate of sBHB mergers as a function of SNR is shown in Figure 5
for 4-year (purple) and 10-year (teal) LISA mission for model D and model DK.

During a 4-year observational run, our results for model D indicate that LISA
and LIGO will detect $27$ sBHBs ${\rm yr^{-1}}$ with SNR $\geq 4$. However,
this number drops to $2$ ${\rm yr^{-1}}$ when we increase the SNR threshold to
10 for model D. Model DK also shows the same pattern but is six times smaller
than model D, i.e $5$ ${\rm yr^{-1}}$ for SNR $\geq 4$ to 0.3 ${\rm yr^{-1}}$
for an SNR $\geq 10$.

If we extend the observational run to 10 years, the detectable multiband rate
increases to $132$${\rm yr^{-1}}$ with an SNR threshold of 4 to $6$ ${\rm
yr^{-1}}$ for SNR $\geq 10$. For model DK we see only $22$ ${\rm yr^{-1}}$ for
SNR $\geq 4$ to $1$ ${\rm yr^{-1}}$ for SNR $\geq 10$.

###  3.2 Expected Number of Multiband sBHBs

Converting the multiband rate from equation 4 to the number of detectable
sBHBs, we used:

$N=\frac{dN}{dzdt}T_{\rm obs},$ (5)

where $\frac{dN}{dzdt}$ is the rate from equation 4 and $T_{\rm obs}$ is the
LISA mission lifetime. Figure 6 displays the number of sBHB that can be
detected by both LIGO and LISA given a 4- and 10-year LISA mission lifetime.
We see once again that the detectable multiband BH rate decreases as the SNR
thresholds increase for both a 4-year and 10-year mission. For a 4-year LISA
run, model D (left panel) exhibits a detectable number of multiband sBHB
mergers that ranges from 6–106 as the SNR threshold decreases from 10 to 4,
while model DK (right panel) yields between 1–18. For a 10-year mission, the
number of sBHB mergers in model D ranges from 59–1321, and from 10–220 for
model DK.

From our results, we conclude that the multiband detection rates (per year)
for a 10-year LISA mission increases by ${\sim}4$x for both the models across
all SNR thresholds when compared to a 4-year LISA mission. This is consistent
with an average increase in SNR of sBHB by around $\sqrt{10/4}$ between the
two LISA mission lifetimes, thus increasing the detection volume by a factor
of ${\sim}4$. Similarly, the total number of multiband detections in a 10-year
LISA mission is ${\sim}10$x higher at SNR threshold of 10 and ${\sim}12$x
higher at SNR threshold of 4 when compared to 4-year LISA mission. The mild
increase at low SNR thresholds compared to high SNR thresholds is due to an
increase in the intrinsic sBHB merger rates with redshift in the local
universe (Mapelli et al., 2017).

Table 1:  Expected number of multiband sBHBs for a range of LISA mission durations. The terms marginal and expected are defined by SNR $\geq 4$ and SNR $\geq 8$, respectively. For each range, the lower value represents model DK, while the higher value represents model D. | LISA mission
---
lifetime
| Marginal
---
sBHB detections
| Confirmed
---
sBHB detections
4 years | $18-106$ | $2-11$
6 years | $53-320$ | $5-32$
8 years | $118-708$ | $11-68$
10 years | $220-1321$ | $20-112$

###  3.3 Comparison to Previous Work

Previous work has been done to estimate sBHB multiband rates and expected
numbers for LISA (Gerosa et al., 2019; Sesana, 2016, 2017; Kyutoku & Seto,
2016; Zhao et al., 2023). Sesana (2016, 2017) presented a semi-analytic data-
driven approach using LIGO observations from O1, assuming two underlying black
hole mass distributions and intrinsic merger rates: a Salpeter mass function
with a merger rate of $100\,{\rm Gpc^{-3}yr^{-1}}$, and log-uniform mass
function with a rate of $35\,{\rm Gpc^{-3}yr^{-1}}$. These higher merger rates
inferred from O1 naturally lead to a larger number of predicted multiband
events compared to our models.

Kyutoku & Seto (2016) examined the multiband detection of GW150914-analogs
assuming an O1 event rate as well, however this work was primarily focused
comparing the effect of different noise curves for eLISA, a precursor to the
mission architecture adopted by LISA. For a 10 year mission, they predicted
that the number eLISA-LIGO multiband sBHB detections ranged from 8–400,
depending on the mission configuration.

Gerosa et al. (2019) updated the data-driven approach of the Sesana work using
LIGO/Virgo/KAGRA observations from O1 and O2; the Salpeter mass function
yielded a merger rate of $57^{+40}_{-25}\,{\rm Gpc^{-3}yr^{-1}}$, while the
log-uniform yielded a rate of $19^{+13}_{-8.2}\,{\rm Gpc^{-3}yr^{-1}}$. For a
10-year LISA mission, their finding of 4-22 multiband LISA+LVK detections
above an SNR of 8 is most consistent with our model DK, while model D is
roughly 6 times higher.

Additionally, Gerosa et al. (2019) created binary population synthesis
simulations with spinning black hole binaries for isolated systems that
evolved through the common-envelope phase. They perform multiple models, each
varying the natal kicks distributed using a Maxwellian distribution with a
dispersion between 0 and 265 km/s, and noted that the expected multiband
detection rates can vary by more than an order of magnitude depending on the
natal kick choice, consistent with our findings.

Zhao et al. (2023) created a Monte Carlo-sampled database of mock sBHB based
on GWTC-3 observations with a rate of $19^{+8.4}_{-8.5}$
$\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$. For a 4-year LISA mission and an SNR
threshold of 8 they predict 3-15 events, which falls within the range of both
of our models.

##  4 Conclusion

In this paper, we aimed to predict the expected rate and signature of stellar
origin binary black holes that will be detectable by both LISA and LVK, with
the hope that the early warning of a merger will help constrain sky location
and help facilitate an alert of an impending merger for electromagnetic and
ground-based GW detectors. Combining two binary population synthesis models
with a cosmological hydrodynamic simulation, we found the multiband rate and
expected number of sBHB detectable as a function of SNR and LISA mission
duration. Table 1 summarizes our findings for four different mission
durations. In this table, we define marginal detections as those with SNR
$\geq 4$, while confirmed sBHB detections must have SNR $\geq 8$. Model DK
produces the minimum value in each column while model D yields the maximum
value. We note that model DK, the binary population synthesis model with
larger supernova natal kicks, yields overall sBHB merger rates that are more
consistent with LVK observations through O3. Although it may be tempting to
conclude that small black hole natal kicks can be ruled out by LVK
observations, the number of competing open questions in binary stellar
evolution – from the strength of core overshooting, to the impact of unstable
mass transfer, to the core-collapse explosion process itself – preclude the
ability to do so. Advances in binary stellar population modeling will allow a
better mapping of the unknown parameter space in binary stellar evolution and
will therefore better enable comparisons to gravitational wave observations
(Spera et al., 2019).

Our work provides a novel approach for building LISA mock data catalogs of
sBHB directly using cosmological simulations. Such synthetic catalogs provide
useful benchmark for LISA data-analysis pipelines (Babak et al., 2010; Baghi,
2022; Lackeos et al., 2023). To maximize the LISA SNR of a given sBHB, we
deliberately modeled those that would merge over a timescale that is roughly
the LISA mission duration; there are certainly sBHBs that will merge in a
ground-based gravitational wave band hundreds of years in the future, but
these sources are much fainter in the LISA band when so far from merger, and
are therefore not ideal multiband sources. We also neglected sBHBs that did
not merge by the end of the simulation; this population of monochromatic, or
nearly monochromatic, sBHBs may well be buried in the confusion foreground of
Galactic compact object binaries (Colpi et al., 2024).

This study considers only the field sBHB formation channel, yet LVK
observations may hint that multiple formation channels could be at play
(Bouffanais et al., 2021). Cosmological simulations of the scale needed to
make universal predictions are as yet unable to resolve AGN disks and dense
stellar systems such as nuclear star clusters and globular clusters. We will
explore building a framework to seed a simulation with sBHB from multiple
formation channels in a future work.

## Acknowledgments

KHB acknowledges and appreciates support from NSF NRT-2125764. M.M.
acknowledges financial support from the European Research Council for the ERC
Consolidator grant DEMOBLACK, under contract No. 770017 (PI: Mapelli) and from
the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC
2181 - 390900948) STRUCTURES. This work used the resources provided by the
Vanderbilt Advanced Computing Center for Research and Education (ACCRE), a
collaboratory operated by and for Vanderbilt faculty at Vanderbilt University.

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102, doi: 10.1103/PhysRevLett.116.061102
  * Abbott et al. (2016a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, ApJ, 818, L22, doi: 10.3847/2041-8205/818/2/L22
  * Abbott et al. (2016b) —. 2016b, Phys. Rev. Lett., 116, 061102, doi: 10.1103/PhysRevLett.116.061102
  * Abbott et al. (2019) —. 2019, Physical Review X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, Physical Review X, 11, 021053, doi: 10.1103/PhysRevX.11.021053
  * Abbott et al. (2023a) Abbott, R., Abbott, T. D., Acernese, F., et al. 2023a, Physical Review X, 13, 011048, doi: 10.1103/PhysRevX.13.011048
  * Abbott et al. (2023b) —. 2023b, Physical Review X, 13, 041039, doi: 10.1103/PhysRevX.13.041039
  * Abbott et al. (2024) —. 2024, Phys. Rev. D, 109, 022001, doi: 10.1103/PhysRevD.109.022001
  * Ajith et al. (2007) Ajith, P., Babak, S., Chen, Y., et al. 2007, Classical and Quantum Gravity, 24, S689, doi: 10.1088/0264-9381/24/19/S31
  * Amaro-Seoane et al. (2017a) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017a, arXiv e-prints, arXiv:1702.00786.  https://arxiv.org/abs/1702.00786
  * Amaro-Seoane et al. (2017b) —. 2017b, arXiv e-prints, arXiv:1702.00786.  https://arxiv.org/abs/1702.00786
  * Askar et al. (2017) Askar, A., Szkudlarek, M., Gondek-Rosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36, doi: 10.1093/mnrasl/slw177
  * Babak et al. (2010) Babak, S., Baker, J. G., Benacquista, M. J., et al. 2010, Classical and Quantum Gravity, 27, 084009, doi: 10.1088/0264-9381/27/8/084009
  * Babak et al. (2023) Babak, S., Caprini, C., Figueroa, D. G., et al. 2023, J. Cosmology Astropart. Phys, 2023, 034, doi: 10.1088/1475-7516/2023/08/034
  * Baghi (2022) Baghi, Q. 2022, arXiv e-prints, arXiv:2204.12142, doi: 10.48550/arXiv.2204.12142
  * Baker et al. (2019) Baker, J., Bellovary, J., Bender, P. L., et al. 2019, arXiv e-prints, arXiv:1907.06482.  https://arxiv.org/abs/1907.06482
  * Banerjee (2017) Banerjee, S. 2017, MNRAS, 467, 524, doi: 10.1093/mnras/stw3392
  * Banerjee et al. (2010) Banerjee, S., Baumgardt, H., & Kroupa, P. 2010, MNRAS, 402, 371, doi: 10.1111/j.1365-2966.2009.15880.x
  * Bartos et al. (2017) Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, ApJ, 835, 165, doi: 10.3847/1538-4357/835/2/165
  * Bassan (2014) Bassan, M. 2014, Advanced interferometers and the search for gravitational waves: Lectures from the first VESF school on advanced detectors for Gravitational Waves (Springer)
  * Belczynski et al. (2016) Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016, Nature, 534, 512, doi: 10.1038/nature18322
  * Belczynski et al. (2002) Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407, doi: 10.1086/340304
  * Bethe & Brown (1998) Bethe, H. A., & Brown, G. E. 1998, ApJ, 506, 780, doi: 10.1086/306265
  * Bouffanais et al. (2021) Bouffanais, Y., Mapelli, M., Santoliquido, F., et al. 2021, MNRAS, 507, 5224, doi: 10.1093/mnras/stab2438
  * Breivik et al. (2016) Breivik, K., Rodriguez, C. L., Larson, S. L., Kalogera, V., & Rasio, F. A. 2016, ApJ, 830, L18, doi: 10.3847/2041-8205/830/1/L18
  * Chan et al. (2018) Chan, M. L., Messenger, C., Heng, I. S., & Hendry, M. 2018, Phys. Rev. D, 97, 123014, doi: 10.1103/PhysRevD.97.123014
  * Chen et al. (2022) Chen, N., Ni, Y., Holgado, A. M., et al. 2022, MNRAS, 514, 2220, doi: 10.1093/mnras/stac1432
  * Colpi et al. (2024) Colpi, M., Danzmann, K., Hewitson, M., et al. 2024, LISA Definition Study Report.  https://arxiv.org/abs/2402.07571
  * Croton (2013) Croton, D. J. 2013, PASA, 30, e052, doi: 10.1017/pasa.2013.31
  * Cutler et al. (2019) Cutler, C., Berti, E., Holley-Bockelmann, K., et al. 2019, Bulletin of the American Astronomical Society, 51, 109, doi: 10.48550/arXiv.1903.04069
  * de Mink & Mandel (2016) de Mink, S. E., & Mandel, I. 2016, MNRAS, 460, 3545, doi: 10.1093/mnras/stw1219
  * Digman & Cornish (2023) Digman, M. C., & Cornish, N. J. 2023, Physical Review D, 108, doi: 10.1103/physrevd.108.023022
  * Ford & McKernan (2022) Ford, K. E. S., & McKernan, B. 2022, MNRAS, 517, 5827, doi: 10.1093/mnras/stac2861
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91, doi: 10.1088/0004-637X/749/1/91
  * Gerosa et al. (2019) Gerosa, D., Ma, S., Wong, K. W. K., et al. 2019, Phys. Rev. D, 99, 103004, doi: 10.1103/PhysRevD.99.103004
  * Giacobbo et al. (2018) Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959, doi: 10.1093/mnras/stx2933
  * Hinshaw et al. (2013) Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19, doi: 10.1088/0067-0049/208/2/19
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974, doi: 10.1111/j.1365-2966.2005.09087.x
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543, doi: 10.1046/j.1365-8711.2000.03426.x
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897, doi: 10.1046/j.1365-8711.2002.05038.x
  * Jani et al. (2020) Jani, K., Shoemaker, D., & Cutler, C. 2020, Nature Astronomy, 4, 260, doi: 10.1038/s41550-019-0932-7
  * Korol et al. (2017) Korol, V., Rossi, E. M., Groot, P. J., et al. 2017, Monthly Notices of the Royal Astronomical Society, 470, 1894, doi: 10.1093/mnras/stx1285
  * Kruckow et al. (2018) Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, MNRAS, 481, 1908, doi: 10.1093/mnras/sty2190
  * Kyutoku & Seto (2016) Kyutoku, K., & Seto, N. 2016, MNRAS, 462, 2177, doi: 10.1093/mnras/stw1767
  * Lackeos et al. (2023) Lackeos, K., Littenberg, T. B., Cornish, N. J., & Thorpe, J. I. 2023, Astronomy & Astrophysics, 678, A123, doi: 10.1051/0004-6361/202347222
  * Lamberts et al. (2019) Lamberts, A., Blunt, S., Littenberg, T. B., et al. 2019, MNRAS, 490, 5888, doi: 10.1093/mnras/stz2834
  * Madau & Fragos (2017) Madau, P., & Fragos, T. 2017, ApJ, 840, 39, doi: 10.3847/1538-4357/aa6af9
  * Mapelli (2016) Mapelli, M. 2016, MNRAS, 459, 3432, doi: 10.1093/mnras/stw869
  * Mapelli et al. (2017) Mapelli, M., Giacobbo, N., Ripamonti, E., & Spera, M. 2017, MNRAS, 472, 2422, doi: 10.1093/mnras/stx2123
  * Marchant et al. (2016) Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50, doi: 10.1051/0004-6361/201628133
  * McKernan et al. (2018) McKernan, B., Ford, K. E. S., Bellovary, J., et al. 2018, ApJ, 866, 66, doi: 10.3847/1538-4357/aadae5
  * Mennekens & Vanbeveren (2014) Mennekens, N., & Vanbeveren, D. 2014, A&A, 564, A134, doi: 10.1051/0004-6361/201322198
  * Portegies Zwart & McMillan (2000) Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17, doi: 10.1086/312422
  * Randall & Xianyu (2021) Randall, L., & Xianyu, Z.-Z. 2021, The Astrophysical Journal, 914, 75, doi: 10.3847/1538-4357/abfd2c
  * Ranjan et al. (2024) Ranjan, S., Jani, K., Nitz, A. H., Holley-Bockelmann, K., & Cutler, C. 2024, arXiv e-prints, arXiv:2406.11926, doi: 10.48550/arXiv.2406.11926
  * Robson et al. (2019) Robson, T., Cornish, N. J., & Liu, C. 2019, Classical and Quantum Gravity, 36, 105011, doi: 10.1088/1361-6382/ab1101
  * Rodriguez et al. (2016) Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016, Phys. Rev. D, 93, 084029, doi: 10.1103/PhysRevD.93.084029
  * Rodriguez et al. (2021) Rodriguez, C. L., Kremer, K., Chatterjee, S., et al. 2021, Research Notes of the American Astronomical Society, 5, 19, doi: 10.3847/2515-5172/abdf54
  * Rodriguez et al. (2015) Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, Phys. Rev. Lett., 115, 051101, doi: 10.1103/PhysRevLett.115.051101
  * Sesana (2016) Sesana, A. 2016, Phys. Rev. Lett., 116, 231102, doi: 10.1103/PhysRevLett.116.231102
  * Sesana (2017) Sesana, A. 2017, in Journal of Physics Conference Series, Vol. 840, Journal of Physics Conference Series (IOP), 012018, doi: 10.1088/1742-6596/840/1/012018
  * Spera et al. (2015) Spera, M., Mapelli, M., & Bressan, A. 2015, MNRAS, 451, 4086, doi: 10.1093/mnras/stv1161
  * Spera et al. (2019) Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889, doi: 10.1093/mnras/stz359
  * Springel (2010) Springel, V. 2010, MNRAS, 401, 791, doi: 10.1111/j.1365-2966.2009.15715.x
  * Stone et al. (2017) Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946, doi: 10.1093/mnras/stw2260
  * Toubiana et al. (2021) Toubiana, A., Sberna, L., Caputo, A., et al. 2021, Physical Review Letters, 126, doi: 10.1103/physrevlett.126.101105
  * van den Heuvel (2019) van den Heuvel, E. P. J. 2019, IAU Symposium, 346, 1, doi: 10.1017/S1743921319001315
  * Vitale (2016) Vitale, S. 2016, Physical Review Letters, 117, doi: 10.1103/physrevlett.117.051102
  * Vogelsberger et al. (2014a) Vogelsberger, M., Genel, S., Springel, V., et al. 2014a, MNRAS, 444, 1518, doi: 10.1093/mnras/stu1536
  * Vogelsberger et al. (2014b) —. 2014b, Nature, 509, 177, doi: 10.1038/nature13316
  * Zhao et al. (2023) Zhao, Y., Lu, Y., Yan, C., Chen, Z., & Ni, W.-T. 2023, MNRAS, 522, 2951, doi: 10.1093/mnras/stad1107
  * Ziosi et al. (2014) Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703, doi: 10.1093/mnras/stu824


Hamiltonian singular value transformation and inverse block encoding

Seth Lloyd,∗1,2 Bobak T. Kiani,1,3 David R.M. Arvidsson-Shukur4

Samuel Bosch,1,3 Giacomo De Palma,5, William M. Kaminsky,1

Zi-Wen Liu,7 Milad Marvian8

1\. Turing, Inc. 2\. Department of Mechanical Engineering, MIT.

3\. Department of Electrical Engineering and Computer Science, MIT. 4. Hitachi
Laboratory, Cambridge UK.

5\. Scuola Normale Superiore, Pisa IT. 6. Department of Physics, MIT. 7.
Perimeter Institute CA.

8\. Department of Physics, Department of Electrical and Computer Engineering,
UNM.

∗ to whom correspondence should be addressed<EMAIL_ADDRESS>

Abstract: The quantum singular value transformation is a powerful quantum
algorithm that allows one to apply a polynomial transformation to the singular
values of a matrix that is embedded as a block of a unitary transformation.
This paper shows how to perform the quantum singular value transformation for
a matrix that can be embedded as a block of a Hamiltonian. The transformation
can be implemented in a purely Hamiltonian context by the alternating
application of Hamiltonians for chosen intervals: it is an example of the
Quantum Alternating Operator Ansatz (generalized QAOA). We also show how to
use the Hamiltonian quantum singular value transformation to perform inverse
block encoding to implement a unitary of which a given Hamiltonian is a block.
Inverse block encoding leads to novel procedures for matrix multiplication and
for solving differential equations on quantum information processors in a
purely Hamiltonian fashion.

The quantum singular value transformation (QSVT) [1] shows that if one can
embed a matrix $A$ in a block of a performable unitary transformation, one can
also implement a transformation with the same singular vectors as $A$, but
whose singular values are polynomial functions of the singular values of $A$.
The quantum singular value transformation is a powerful quantum algorithm that
encompasses other quantum algorithms such as quantum simulation and matrix
inversion as subcases. The methods of the QSVT are based on the technique of
qubitization [2], a state-of-the-art method for simulating a Hamiltonian
evolution given the ability to perform block encoding.

This paper presents a purely Hamiltonian version of the quantum singular value
transformation. We show that if one has the ability to apply a Hamiltonian $H$
of which $A$ is a block, then we can perform the full quantum singular value
transformation using the Quantum Alternating Operator Ansatz (generalized
QAOA): we alternate application of the Hamiltonian $H$ with a second, readily
applied Hamiltonian, for a chosen set of times. Our motivation in developing
Hamiltonian QSVT is simple: quantum information processing is generically
implemented via the application of semiclassical control fields, which apply
time-dependent Hamiltonians to the quantum information bearing degrees of
freedom. In the gate model of quantum computing, these control fields are used
to implement unitary quantum logic gates, which are applied to the
qubits/qudits/qumodes in the quantum computer. But many more possible
Hamiltonians can be applied, including global Hamiltonians, as in the Quantum
Approximate Optimization Algorithm (original QAOA) [3]. In general, the
ability to apply a set of Hamiltonians governed by the application of time
varying control fields allows one to apply any effective Hamiltonian in the
algebra generated by the set. Different coherent quantum information
processors, e.g., coherent diabatic quantum annealers [4], allow the
application of the Hamiltonian quantum singular value transformation to the
set of Hamiltonians that are ‘native’ to those processors – i.e., those that
can be efficiently applied. By introducing a purely Hamiltonian version of the
quantum singular value transformation that bypasses the gate model of quantum
computing, we hope to expand the power of near term quantum information
processors.

Main result:

In this paper we show that if one can apply a Hamiltonian of which a matrix
$A$ is an off-diagonal block, we can deterministically implement the unitary
transformation

$U_{f}=i\pmatrix{\sqrt{I-f(A^{\dagger})f(A)}&f(A^{\dagger})\cr
f(A)&-\sqrt{I-f(A)f(A^{\dagger})}\cr},$ $None$

where $f(A)$ is the matrix with the same singular vectors as $A$, and singular
values that are a function $f$ of the singular values of $A$.

$U_{f}$ is easily verified to be unitary as long as $f(A^{\dagger})f(A)\leq
I$. (For example, expand the square root as a power series.) Note that $U$ is
anti-Hermitian as well as unitary, and so has eigenvalues $\pm i$. Like the
unitary QSVT, the Hamiltonian singular value transformation takes time
$O(\sigma_{min}^{-1}\log(1/\epsilon))$ to perform the transformation to
accuracy $\epsilon$, where $\sigma_{min}$ is the smallest singular value of
$A$.

Setting $f(x)=x$ yields inverse block encoding, or inverse qubitization [2]:
the ability to apply a Hamiltonian of which a matrix is a block allows one to
apply a unitary of which the matrix is a block. Below, we apply inverse block
encoding to give Hamiltonian versions of matrix multiplication and the quantum
solution of differential equations.

Preliminaries:

Assume the ability to apply a Hamiltonian $\pm H$, where

$H=\pmatrix{\bullet&A^{\dagger}\cr A&\bullet\cr},$ $None$

and the on-diagonal blocks are arbitrary. We assume that $A^{\dagger}A\leq I$.
We also assume the ability to apply the Hamiltonian $\pm Z$, where

$Z\equiv\pmatrix{I&0\cr 0&-I\cr}.$ $None$

As shown in [5], by applying $\pm H$ in alternation with $\pm Z$, we can
average out/refocus the on-diagonal terms of $H$ in time $O(1)$, so without
loss of generality we take $H$ to be of the form

$H=\pmatrix{0&A^{\dagger}\cr A&0\cr}.$ $None$

Write $A=\sum_{j}\sigma_{j}|\ell_{j}\rangle\langle r_{j}|$, where $\sigma_{j}$
are the singular values of $A$ and $|r_{j}\rangle$, $|\ell_{j}\rangle$ are the
corresponding right and left singular vectors. The Hamiltonian $H$ acts on the
Hilbert space ${\cal H}_{R}\oplus{\cal H}_{L}$ consisting of the direct sum of
the the Hilbert spaces for the right and left singular vectors. The
eigenvectors of $H$ take the form

$H\pmatrix{|r_{j}\rangle\cr\pm|\ell_{j}\rangle\cr}=\pm\sigma_{j}\pmatrix{|r_{j}\rangle\cr\pm|\ell_{j}\rangle\cr}.$
$None$

That is, within the two-dimensional subspace ${\cal H}_{j}$ spanned by
$\pmatrix{|r_{j}\rangle\cr\pm|\ell_{j}\rangle\cr}$, $H$ acts as $\sigma_{j}X$,
where $X$ is the $x$ Pauli matrix in this subspace.

Note that we can also apply the transformation

$e^{i\phi Z/2}e^{-iHt}e^{-i\phi Z/2}\equiv e^{-iG_{\phi}t},$ $None$

where

$G_{\phi}\equiv\pmatrix{0&e^{i\phi}A^{\dagger}\cr e^{-i\phi}A&0\cr}.$ $None$

$G_{\phi}$ acts within ${\cal H}_{j}$ as $\sigma_{j}$ times a Pauli matrix for
a rotation about an axis at angle $\phi$ in the $xy$ plane.

Hamiltonian quantum singular value transformation:

From the results of the previous section, we can apply transformations of the
form

$e^{-iG_{\phi_{k}}t_{k}}\ldots e^{-iG_{\phi_{1}}t_{1}}.$ $None$

Equation (8) shows that we can generate unitary transformations that act in
each of the ${\cal H}_{j}$ as polynomials in $\cos(\sigma_{j}t_{i})$,
$\sin(\sigma_{j}t_{i})$. The coefficents of these polynomials are determined
by tuning the $\phi_{i}$. Qubitization [2] shows how to generate a broad class
of such polynomial unitary transformations. Note that because of the
Hamiltonian nature of our construction, we have greater flexibility than the
transformations applied in qubitization [2] and the quantum singular value
transformation [1]: there, because $H$ is applied as a block of a unitary, the
times $t_{i}$ are all fixed to be equal. Here, we can vary the times $t_{i}$
as well as the angles $\phi_{i}$. This feature could prove useful for
optimizing the transformation using variational methods.

The theorems of [1-2] show that the ability to perform transformations of the
form (8) allows us to apply transformations of the form

$U=\bigoplus_{j}\pmatrix{P(\cos\sigma_{j})&i\sin\sigma_{j}Q(\cos\sigma_{j})\cr
i\sin\sigma_{j}Q^{*}(\cos\sigma_{j})&P^{*}(\cos\sigma_{j})},$ $None$

where the degree of the polynomial $P$ is less than or equal to $k$, and the
degree of the polynomial $Q$ is less than or equal to $k-1$. The polynomials
$P$ and $Q$ can be determined at will by choosing the angles $\phi_{i}$,
subject to the constraint that $P$ has parity $k$ mod 2, $Q$ has parity $k-1$
mod 2, and that $|P(\cos\sigma)|^{2}+\sin^{2}\sigma|Q(\cos\sigma)|^{2}=1$. The
first of these constraints comes from the product form of equation (8), and
the second from unitarity.

To perform the Hamiltonian singular value transformation, we choose $Q$ so
that

$\sin\sigma_{j}Q^{*}(\cos\sigma_{j})=f(\sigma_{j})\pm\epsilon$ $None$

Equivalently, $\sqrt{1-x^{2}}Q^{*}(x)\approx f(\arccos(x))$: as in [1-2] we
are effectively decomposing the function $f(\arccos(x))$ in terms of Chebyshev
polynomials. For suitably smooth functions $f$ the decomposition achieves
accuracy $\epsilon$ for Chebyshev polyomials of degree $O(\log(1/\epsilon))$.
Similarly, take

$-iP(\cos(\sigma_{j}))=\sqrt{I-f^{*}(\sigma_{j})f(\sigma_{j})}\pm\epsilon.$
$None$

To within accuracy $\epsilon$, the resulting unitary transformation is

$U_{f}=i\pmatrix{\sqrt{I-f(A^{\dagger})f(A)}&f(A^{\dagger})\cr
f(A)&-\sqrt{I-f(A)f(A^{\dagger})}\cr},$ $None$

which is the desired Hamiltonian singular value transformation of equation
(1).

Inverse block encoding is obtained by taking $f(\sigma)=\sigma$, which gives

$U=i\pmatrix{\sqrt{I-A^{\dagger}A}&A^{\dagger}\cr
A&-\sqrt{I-AA^{\dagger}}\cr}.$ $None$

Because the derivative of $\arccos$ diverges at $\pm 1$, we require
$|\sigma|<1-\delta$ in order to obtain an error scaling of
$\epsilon=O(e^{-\sqrt{2\delta}k})$ in the $k$-th order Chebyshev polynomial
expansion of $\arccos(x)$ [6-7].

The error scaling for the Hamiltonian quantum singular value transformation is
the same as for the unitary QSVT, and the number of alternating steps required
to perform the transformation goes as
$\\#=O(\sigma_{min}^{-1}\log(1/\epsilon))$. In comparison with performing the
same procedure using the generalized quantum linear systems algorithm [5,8],
which relies on quantum phase estimation, we see that the dependence of the
time has gone from $1/\epsilon$ to $\log(1/\epsilon)$.

Since in practice the alternating Hamiltonians are applied by turning on and
off semiclassical control fields, we also need to address the question of how
accurately these control fields need to be applied. If each Hamiltonian
application consists of a time-dependent term of the form $\gamma(t)H$, then
the key figure of merit is the total integral of $\gamma(t)$ over the period
during which the control is turned on and off: this integral needs to be
precise to accuracy $O(\epsilon/\\#)$ in order for the accumulated error over
the overall alternating sequence to be less than $\epsilon$. Note that in
existing high-coherence quantum information processors such as ion traps, the
semiclassical control error is the dominant source of error over the quantum
computation. Current superconducting and ion-trap quantum computers have pulse
control errors on the order of $2*10^{-3}-5*10^{-3}$, with a next-generation
goal of $1*10^{-3}$. Given the favorable error scaling, the Hamiltonian QSVT
should be performable on current NISQ devices.

The Hamiltonian quantum singular value transformation is implemented solely
through the alternating application of the Hamiltonian $H$ and the Hamiltonian
$Z$. That is, it is an application of the Quantum Alternating Operator Ansatz
or generalized QAOA. As in the original QAOA (Quantum Approximate Optimization
Algorithm) [3], this feature makes the Hamiltonian QSVT ripe for the
application of variational methods: we set some task for the unitary $U_{f}$
to accomplish and try to find an optimal $f$ by varying the parameters
$t_{i},\phi_{i}$ in equation (8).

Applications:

Inverse block encoding provides a novel method for quantum matrix
multiplication [9-10]: from equation (13) we have

$U\pmatrix{|\psi\rangle\cr
0\cr}=i\pmatrix{(\sqrt{I-A^{\dagger}A})|\psi\rangle\cr A|\psi\rangle\cr},$
$None$

which allows one to implement $A|\psi\rangle$ with probability

$\langle\psi|A^{\dagger}A|\psi\rangle.$ $None$

Amplitude amplification provides a quadratic improvement, allowing the state
$A|\psi\rangle$ to be prepared in time $\sqrt{\psi|A^{\dagger}A|\psi\rangle}$.

When $A$ is square, applying this procedure first to the first and second
components of the state $(|\psi\rangle,0,\ldots,0)^{T}$, then to the second
and third components, then to the third and fourth, etc. allows one
deterministically to construct the state

$\pmatrix{(\sqrt{I-A^{\dagger}A})|\psi\rangle\cr(\sqrt{I-A^{\dagger}A})A|\psi\rangle\cr\ldots\cr(\sqrt{I-A^{\dagger}A})A^{n-1}|\psi\rangle\cr
A^{n}|\psi\rangle\cr}=\sum_{k=0}^{n-1}(\sqrt{I-A^{\dagger}A})A^{k}|\psi\rangle|k\rangle+A^{n}|\psi\rangle|n\rangle.$
$None$

Measuring the value of the second register yields the state
$A^{n}|\psi\rangle$ with probability

$\langle\psi|{A^{\dagger}}^{n}A^{n}|\psi\rangle.$ $None$

Again, amplitude amplification provides a quadratic speed up in the
preparation of this state.

Note that the requirement $A^{\dagger}A\leq I$ means that the probability of
obtaining $A^{n}|\psi\rangle$ will be asymptotically exponentially suppressed
except for components corresponding to the singular value $\sigma_{j}=1$.
Since $\sigma_{j}=1$ corresponds to subspaces on which $H$ acts in a unitary
fashion, this method allows us to identify ‘hidden’ unitaries in $H$.

When $A=I+B\Delta t$, then the state (16) yields the Euler forward solution of
the differential equation

${d|\psi\rangle\over dt}=B|\psi\rangle.$ $None$

The requirement that

$A^{\dagger}A=I+(B+B^{\dagger})\Delta t+B^{\dagger}B\Delta t^{2}\leq I$ $None$

means that the eigenvalues of $B$ must have negative real part, so that the
differential equation (18) is dissipative. If $B$ is sparse or low rank, then
we can implement the Hamiltonian $H$ via conventional methods of quantum
Hamiltonian simulation [11-13]. This method of solving the forward
differential equation allows one to obtain the (unnormalized) state
$|\psi(t)\rangle=e^{Bt}|\psi(0)\rangle$ that is the solution to the
differential equation (18) by projecting onto the final component of the state
in (16), which succeeds with probability $\langle\psi(t)|\psi(t)\rangle$.

That is, one can use inverse block encoding to solve dissipative differential
equations on a quantum computer without having to resort to matrix inversion,
a potential savings in computational complexity over conventional quantum
differential equation solvers [14-15].

History states:

The Hamiltonian version of the quantum singular value transformation can be
used to perform quantum matrix inversion to apply
$(\sqrt{I-A^{\dagger}A})^{-1}$ to the first $n$ components of the vector (16)
to obtain the state

$\pmatrix{|\psi\rangle\cr A|\psi\rangle\cr\ldots\cr
A^{n}|\psi\rangle}=\sum_{k=0}^{n}A^{k}|\psi\rangle|k\rangle.$ $None$

In the case that $A=I+B\Delta t$, this state is the ‘history state’ for the
solution of the differential equation (18). The matrix inversion is non-
deterministic, and succeeds with probability $O(1/\tilde{\kappa})$, where
$\tilde{\kappa}$ is the condition number of $\sqrt{I-A^{\dagger}A}$.

Multiple Hamiltonians:

When building quantum computers, we implement unitary transformations by
applying time dependent Hamiltonians $\sum_{j}\gamma_{j}(t)H_{j}$ governed by
semiclassical control fields $\gamma_{j}(t)$. This ability allows us to
implement effective Hamiltonians in the Lie algebra ${\cal A}$ generated by
the $H_{j}$, and to implement the corresponding embedding unitaries. The
Hamiltonian QSVT then allows us deterministically to implement the $U_{f}$ of
equation (1) for any $A$ that can be embedded in a block of a Hamiltonian
$\in{\cal A}$. Note that even systems where the applied Hamiltonians are quite
simple, as in QAOA [3] where one alternates the application of an Ising
Hamiltonian with identical rotations of the qubits about the $X$ axis, the
algebra generated by the Hamiltonians is typically the full algebra of all
possible Hamiltonians – QAOA is universal for quantum computation [16]. The
set of effective Hamiltonians that can be generated efficiently by any given
coherently controllable quantum system depends on the nature of the control
Hamiltonians. Since Hamiltonian QSVT can be applied to any coherently
controllable Hamiltonian system, we hope that our results will encourage
experimentalists to investigate novel architectures for Hamiltonian quantum
information processing. For example, Hamiltonian QSVT is a natural procedure
to implement on the next generation of coherent diabatic quantum annealers.

Discussion:

This paper showed how to perform a Hamiltonian version of the quantum singular
value transformation solely by the alternating application of Hamiltonians for
different lengths of time. The method is an example of the quantum alternating
operator ansatz (generalized QAOA), and allows the direct construction of the
unitary transformations of the QSVT, without the application of quantum logic
gates. Our goal in presenting this construction is straightforward: quantum
information processing is typically implemented by the application of
semiclassical control fields to implement time-varying Hamiltonians, which can
be applied globally as well as locally. For experimentally accessible
Hamiltonians, the Hamiltonian QSVT represents a direct and simple way to
perform the powerful quantum singular value transformation.

As a simple example of Hamiltonian QSVT, we showed how to perform inverse
block encoding: given the ability to apply a Hamiltonian, deterministically
implement a unitary transformation in which the Hamiltonian is embedded as the
off-diagonal blocks. Inverse block encoding provides novel methods for quantum
matrix multiplication and the quantum solution of forward differential
equations. Given the ability to apply multiple Hamiltonians, inverse block
encoding allows one to embed any transformation in the algebra generated by
the set of accessible Hamiltonians in a deterministically constructed unitary
transformation.

References

[1] A. Gilyén, Y. Su, G.H. Low, N. Wiebe, STOC 2019: Proceedings of the 51st
Annual ACM SIGACT Symposium on Theory of Computing, 193–204 (2019).

[2] G.H. Low, I.L. Chuang, Quantum 3, 163 (2019); arXiv: 1610.06546v3.

[3] E. Farhi, J. Goldstone, S. Gutmann, ‘A quantum approximate optimization
algorithm,’ arXiv: 1411.4028 (2014).

[4] E.J. Crosson, D.A. Lidar, ‘Prospects for Quantum Enhancement with Diabatic
Quantum Annealing,’ arXiv: 2008.09913 (2020).

[5] S. Lloyd, S. Bosch, G. De Palma, B. Kiani, Z.-W. Liu, M. Marvian, P.
Rebentrost, D.M. Arvidsson-Shukur ‘Quantum polar decomposition algorithm,’
arXiv: 2006.00841 (2020).

[6] J.P. Boyd, Chebyshev and Fourier Spectral Methods,, Dover Books on
Mathematics, Dover (2001).

[7] A. Gil, J. Segura, N.M. Temme, Numerical Methods for Special Functions,
Society of Industrial and Applied Mathematics (2007).

[8] A.W. Harrow, A. Hassidim, S. Lloyd, Phys. Rev. Lett. 103, 150502 (2009).

[9] C. Shao, ‘Quantum Algorithms to Matrix Multiplication,’ arXiv: 1803.01601.

[10] L. Zhao, Z. Zhao, P. Rebentrost, J. Fitzsimons ‘Compiling basic linear
algebra subroutines for quantum computers,’ arXiv: 1902.10394.

[11] S. Lloyd, Science 273, 1073-1078 (1996).

[12] I. M. Georgescu, S. Ashhab, F. Nori, Rev. Mod. Phys. 86, 153 (2014).

[13] A.M. Childs, A. Ostrander, Y. Su, Quantum 3, 182 (2019); arXiv:
1805.08385v2.

[14] D.W. Berry, J. Phys. A: Math. Theor. 47 105301 (2014).

[15] D.W. Berry, A.M. Childs, A. Ostrander, G. Wang, Comm. Math. Phys. 356,
1057-1081 (2017).

[16] S. Lloyd, ‘Quantum approximate optimization is computationally
universal,’ arXiv: 1812.11075 (2018).