# The Central Kinematics and Black Hole Mass of 4C+37.11

Tirth Surti  Department of Physics/Kavli Institute for Particle Astrophysics
and Cosmology, Stanford University, Stanford, California 94305, USA Roger W.
Romani<EMAIL_ADDRESS>Roger W. Romani Department of Physics/Kavli
Institute for Particle Astrophysics and Cosmology, Stanford University,
Stanford, California 94305, USA  Julia Scharwächter Gemini Observatory/NSF’s
NOIRLab, 670 N’Aohoku Place, Hilo, Hawaii, USA  Alison Peck Visiting Principal
Research Scientist, Department of Astronomy, University of Maryland, 4296
Stadium Dr., College Park, MD 20742-2421  Greg B. Taylor Department of Physics
and Astronomy, University of New Mexico, Albuquerque, NM 87131, USA

###### Abstract

We report on IFU measurements of the host of the radio source 4C+37.11. This
massive elliptical contains the only resolved double compact nucleus at pc-
scale separation, likely a bound supermassive black hole binary (SMBHB).
$i$-band photometry and GMOS-N IFU spectroscopy show that the galaxy has a
large $r_{b}=1.5^{\prime\prime}$ core and that the stellar velocity dispersion
increases inside of a radius of influence $r_{\rm SOI}\approx
1.3^{\prime\prime}$. Jeans Anisotropic Modeling analysis of the core infers a
total SMBHB mass of $2.8^{+0.8}_{-0.8}\times 10^{10}M_{\odot}$, making this
one of the most massive black hole systems known. Our data indicate that there
has been significant scouring of the central kpc of the host galaxy.

Elliptical Galaxies, Core Dynamics, Black Hole Binaries

##  1 Introduction

Maness et al. (2004) found that the radio galaxy 4C+37.11 (0402+379) contains
two central, compact, flat-spectrum, variable components (designated C1 and
C2) with a VLBA-measured separation 7.3 pc at $z=0.055$ and argued that this
is a gravitationally bound SuperMassive Black Hole Binary (SMBHB). Subsequent
VLBA observations by Rodriguez et al. (2006) have bolstered that claim and
even detected a possible relative proper motion of the binary components
(Bansal et al., 2017). The binary’s host is remarkable, being an extremely
massive elliptical embedded in a bright X-ray halo (Romani et al., 2014;
Andrade-Santos et al., 2016) which represents a cluster’s worth of stars and
mass in a relatively isolated galaxy. Such ‘fossil clusters’ are believed to
be the product of multiple major mergers. We report here on a kinematic study
of the central region of this galaxy, which supports a scenario where the
source contains a supermassive black hole pair whose present (and former)
interactions have scoured the galaxy core. The estimated SMBHB mass is very
large, the second largest kinematically measured value in the local universe.

We assume a $H_{0}=69.6$, $\Omega_{M}=0.286$, $\Omega_{\Lambda}=0.714$
cosmology (Bennett et al., 2014), so at the redshift of $z=0.055$, 4C+37.11
has a luminosity distance $D_{L}=246.9$ Mpc, angular diameter distance
$D_{A}=221.9$ Mpc, and a projected scale of $1.076$ kpc/arcsec.

Figure 1: Radial velocity ($V$) and velocity dispersion ($\sigma$) maps of
pPXF measurements in the Voronoi bins. The origin is the AGN location,
determined from the wings of the [NII] emission lines. Solid contours mark
0.9$\times$,0.8$\times$,…0.2$\times$ the peak surface brightness and dashed
lines separate quadrants along the major and minor axes.

Figure 2: Sample spectral fits (red) and residuals for S/N $\sim 20$ Voronoi
bins at three different radii. For plotting purposes, the spectra were
rebinned to $61.2$ km/s, 4$\times$ the width used for the kinematic
extraction. Grey bars show masked pixels from chip gaps and from 2-$\sigma$
iterative clipping.

##  2 Observations, Data Reduction and Measurement

Our kinematic study uses two observations. A 2008 December 12 WIYN 3.6m 2200 s
OPTIC $i$ band image presented in Romani et al. (2014) has image quality FWHM
$<0.6^{\prime\prime}$ and measures the galaxy surface brightness profile. For
spectroscopy we use a 2015 December 1 Gemini 8.1 m $6\times 1200$ s GMOS-N IFU
observation of the galaxy center, downloaded from the Gemini Science archive.
These observations used the B600 grating, covering 6300–9195Å at dispersion
0.47 Å/pixel. Conditions were good, with low airmass throughout. The IFU feeds
a fiber array with a $0.22^{\prime\prime}$ pitch, and with the observations
dithered on the sky, IRAF GMOS pipeline processing calibrated and cleaned the
spectra, assembling a data cube covering $4.5^{\prime\prime}\times
6.2^{\prime\prime}$ at $0.1^{\prime\prime}$ scale, centered on the AGN, with
reduced net exposure at the edges. Measurements of images of the PSF-
distributed line wings of the AGN core find an angular resolution of
$0.64^{\prime\prime}$ (FWHM) and sky lines give a median effective spectral
resolution of $2.0$ Å FWHM.

For our kinematic study, we fit the emission line-free range 7810–9190Å ,
masking a small range near 8230Å affected by a chip gap. We log rebin the
spectra over this range to a velocity scale of $15.3$ km/s. Reliable
measurements of the line-of-sight velocity distribution (LOSVD) require a
signal-to-noise (S/N) of $\sim 20$/resolution element. Our data cube provides
only S/N $\sim$ 7/spaxel in the galaxy center, falling to S/N $\sim 1$ in the
more poorly exposed outskirts of the data cube. Thus we use Voronoi binning
(Cappellari & Copin, 2003) to combine the spaxels. With a target S/N $\sim 20$
at the rebinned velocity scale we get 31 bins (Fig. 1); edge bins generally do
not reach the S/N target. Imperfectly excised cosmic ray events and poor sky
subtraction affect many of the edge pixels with few IFU exposures. We
therefore used iterative sigma clipping to excise outlier pixels from each
wavelength in the Voronoi tiles. The means of the clipped tile wavelength
measurements form the final spectra, with their dispersion as the associated
error vector. Velocity structure is measured with the Penalized PiXel Fitting
method (pPXF) of Cappellari (2023), which convolves template stellar spectra
with a Gauss-Hermite decomposition of the LOSVD. With the limited S/N we fit
only $V$ and $\sigma$, skipping the higher order $h_{3}$ and $h_{4}$ moments.

For the stellar templates, we used a combination of G0–M2 giants from the
Indo-US Stellar Library (Valdes et al., 2004). These templates have a slightly
finer scale of 0.44 Å/pixel, with a median spectral resolution of $\sim 1$ Å.
To match the templates to the host galaxy spectra, we convolve the templates
using the difference in the spectral resolutions, following Section 2.2 in
(Cappellari, 2017). pPXF uses additive and multiplicative polynomials in the
fitting model for both the target and template spectra to reduce sensitivity
to imperfect continuum fluxing and veiling components. We varied the order of
these polynomials, finding that both additive and multiplicative terms
required order $\geq 4$ to achieve stable fits. Here we use 4th-order additive
and 6th-order multiplicative Legendre polynomials (4A+6M) as fiducial. Sample
spectral fits at three different radii are shown in Figure 2.

Radial velocity and velocity dispersion maps from the Voronoi tiles are shown
in Figure 1. We detect little if any coherent rotation, typical for massive
ellipticals. However, given the large size of the Voronoi bins, we cannot
exclude a kinematically distinct core, sometimes seen among non-regular
rotators (Cappellari, 2016). We do see a substantial Keplerian-like peak in
$\sigma(r)$ which flattens out at $\sim$300 km/s beyond $\sim
2^{\prime\prime}$. Although we lack kinematic data at the half-light radius
($R_{e}\sim 12^{\prime\prime}$, see below in the discussion of the Multi-
Gaussian Expansion), we can follow Saglia et al. (2016, App. A) and estimate
the effective velocity dispersion by integrating the surface brightness-
weighted $V_{rms}$ to our outermost bin, finding $\sigma_{e}=332\pm 11$ km/s.
This likely overestimates $\sigma_{e}$, as over a third of our bins are
affected by the central velocity cusp.

For dynamical modeling, we also require the host surface brightness profile.
This was measured from the WIYN $i^{\prime}$ image using the Multi-Gaussian
Expansion (MGE) method of Cappellari (2002). We first masked stars and
galaxies near 4C+37.11. We then obtained the image PSF from 15 unsaturated
stars near the host, using the photutils effective PSF (ePSF) routine (Bradley
et al., 2023). In host fitting we assume axisymmetry, with the surface
brightness modeled as a sum of elliptical Gaussian components aligned along
the direction of the photometric major axis ($x^{\prime}$):

$\Sigma(x^{\prime},y^{\prime})=\sum_{j=1}^{N}S_{j}{\rm
exp}\left[-\frac{1}{2\sigma_{i}}\left(x^{\prime 2}+\frac{y^{\prime
2}}{q_{i}^{2}}\right)\right].$ (1)

Here we have $N$ Gaussian components with surface brightness $S_{j}$ and axial
ratio $q_{j}$. The position angles of the Gaussians were fixed at $PA_{\rm
MGE}=72.4^{\circ}$ and $q_{j}$ had a flat prior 0.75–1. Before fitting we
convolve (1) with a MGE expansion of the WIYN PSF. The MGE fits for the PSF
and the i’ band image are summarized in Tables 1 and 2, respectively. Figure 3
shows the MGE fits for the IFU region and the entire region used in mass
modeling. While somewhat triaxial at large radii, the axisymmetric assumption
is adequate over the IFU kinematics region. Using the MGE, with the routine
mge_half_light_isophote of the package JamPy of (Cappellari, 2008) we find a
half-light radius for 4C+37.11 to be $R_{e}=12.1^{\prime\prime}.$

Figure 3: $i^{\prime}$ image isophote contours (black) and MGE fit ellipses (red). Top: the full region, Bottom: the IFU region with kinematic measurements. The yellow regions correspond to masked pixels in the MGE fits. $j$ | $S_{j}$ |  $\sigma_{j}$ ($\mathrm{arcsec}$) | $q_{j}$
---|---|---|---
1 | 0.23 | 0.15 | 0.71
2 | 1.38 | 0.27 | 0.93
3 | 0.362 | 0.36 | 1.00
4 | 0.10 | 0.56 | 0.96
5 | 0.02 | 0.94 | 1.00
6 | 0.002 | 2.1 | 1.00
Table 1: WIYN PSF MGE $j$ |  $\log(S_{j})$ ($L_{\odot}$ $\mathrm{pc}^{-2}$) |  $\sigma_{j}$ ($\mathrm{arcsec}$) | $q_{j}$
---|---|---|---
1 | 4.18 | 0.11 | 0.75
2 | 3.18 | 1.25 | 0.75
3 | 2.83 | 3.02 | 0.76
4 | 2.29 | 6.38 | 0.75
5 | 2.01 | 17.8 | 0.75
Table 2: PSF-Convolved 4C+37.11 MGE

It was noted in Romani et al. (2014) that the host surface brightness profile
flattened within $\sim 1.5$ kpc from an outer $n=4$ Sérsic profile, indicting
that 4C+37.11 is a cored elliptical. We measure the core radius $r_{b}$ by
fitting a core-Sérsic model, with a point source excess for the AGN light, to
the $i^{\prime}$ band image, using the reconstructed PSF for model convolution
and minimizing with the 2D fitting software imfit (Erwin, 2015). The core-
Sérsic profile is described as

$I(r)=I^{\prime}\left[1+\left(\frac{r_{b}}{r}\right)^{\alpha}\right]^{\frac{\gamma}{\alpha}}\mathrm{exp}\left[-b_{n}\left(\frac{r^{\alpha}+r_{b}^{\alpha}}{r_{e}^{\alpha}}\right)^{\frac{1}{\alpha
n}}\right]$ (2)

where

$I^{\prime}=I_{b}2^{-\gamma/\alpha}\mathrm{exp}\left[b_{n}\left(2^{1/\alpha}r_{b}/r_{e}\right)^{1/n}\right]$
(3)

as per Graham et al. (2003). For our fit we fix the sharpness parameter at
$\alpha=100$ (Graham et al., 2003; Trujillo et al., 2004) and derive the fit
parameters in Table 3, where we have converted intensity at the break radius
to a surface brightness $S_{b,i^{\prime}}$. We do indeed see that the outer
galaxy has $n\approx 4$ and find a core radius $r_{b}=1.60\pm 0.05$ kpc. For
$\alpha>10$, we find that $r_{b}$ varies by $<2\sigma$ from this value.

Core Sérsic Parameters | Value
---|---
PA (deg) | $81.9\pm 0.3$
$\epsilon$ | $0.293\pm 0.003$
$S_{b,i^{\prime}}$ (mag $\mathrm{arcsec}^{-2}$) | $18.26\pm 0.03$
$n$ | $3.92\pm 0.08$
$r_{e}$ (arcsec) | $19.0\pm 0.5$
$r_{b}$ (arcsec) | $1.49\pm 0.05$
$\gamma$ | $0.40\pm 0.03$
Point Source Parameters | Value
$i^{\prime}_{tot}$ (mag) | $20.2\pm 0.3$
Table 3: 2D Core-Sérsic + Point Source model fit parameters. Magnitudes
$S_{b,i^{\prime}}$ and $i^{\prime}$ are corrected for an estimated
$A_{i^{\prime}}=2.2$ mag extinction.

##  3 Dynamical Modeling

Given the low S/N kinematics data extending only to $\sim 3^{\prime\prime}$,
we do not attempt full Schwarzschild orbit-superposition modeling. Instead, we
conducted axisymmetric Jeans Anisotropic Modeling (JAM) using the software
JamPy (Cappellari, 2008), which does not require large-radii kinematics to
achieve good constraints (Simon et al., 2024) and produces BH masses
consistent with Schwarzschild analyses (Cappellari et al., 2010).

Unlike the Schwarzschild method, axisymmetric JAM makes a number of
assumptions about the distribution of stellar velocities. In particular, the
velocity ellipsoid is either assumed to be aligned with spherical
$(r,\phi,\theta)$ or cylindrical $(R,\phi,z)$ coordinates, neither of which
can fully represent a real galaxy. Spherically-aligned velocity ellipsoids are
typically good approximations to slow rotator ellipticals, which are weakly
triaxial and have spherical isophotes in the inner half-light radius (Simon et
al., 2024). Cylindrically-aligned velocity ellipsoids are typically applied to
fast rotator massive ellipticals (Cappellari, 2008) but have also been assumed
for slow and non-regular rotators, as in the ATLAS3d survey (Cappellari et
al., 2013). Here, we compare JAM models using both a cylindrically-aligned and
spherically-aligned velocity ellipsoid. In the cylindrical case, we assume
radial anisotropy $\beta_{z}=1-\sigma_{z}^{2}/\sigma_{R}^{2}$ is constant with
$\sigma_{R}^{2}$ = $b\sigma_{z}^{2}$. In the spherical case, radial anisotropy
$\beta=1-\sigma_{\theta}^{2}/\sigma_{r}^{2}$ is constant. The JAM models
include no tangential anisotropy
$\gamma=1-\sigma_{\phi}^{2}/\sigma_{R}^{2}=0$, since the velocity second
moment predictions are independent of $\gamma$; nonzero tangential anisotropy
is not usually needed to model the overall kinematics of real galaxies
(Cappellari, 2016).

###  3.1 JAM Models

For both the spherically and cylindrically aligned cases, we incorporate a
Self-Consistent Model, assuming that the total mass density follows the
surface brightness profile. In particular, the projected mass density can be
obtained by multiplying each Gaussian component in (1) by a constant mass-to-
light ratio $\Upsilon_{i^{\prime}}$. The parameters for this model are then
$\Upsilon_{i^{\prime}}$, the radial anisotropy parameter
($\sigma_{z}/\sigma_{R}$ for cylindrical or $\sigma_{\theta}/\sigma_{r}$ for
spherical), inclination $i$, and the black hole mass $M_{\bullet}$. Here
inclination $i$ is the angle between the galaxy’s symmetry axis and the line-
of-sight direction and $M_{\bullet}$ represents the total mass of the binary
which is modeled by a Keplerian point mass potential. We represent the PSF of
the kinematic observations with a single circular Gaussian of dispersion
$0.27^{\prime\prime}$ (as described in Section 2). We do not incorporate a
model for the contribution of the dark halo’s potential due to the limited S/N
and radial range of the kinematic data; tests with a NFW halo give no halo
constraints and negligible changes to the other fit parameters.

At S/N $\sim 20$ the central bin has a large
$V_{rms}=\sqrt{V^{2}+\sigma^{2}}\approx 532$ km/s (Fig. 4). Since this has a
strong effect on the black hole mass fit, we tested for sensitivity and
systematic bias by optionally dropping the innermost bin. We also fit using a
lower Voronoi target S/N $=18$, for additional bins across the central
$0.5^{\prime\prime}$. In all cases the $V_{rms}$ rise from the SMBHB is still
apparent.

S/N | Velocity Ellipsoid |  | Central
---
Kin. Bin
$i\text{ (deg)}$ | $\sigma_{\theta}/\sigma_{r}$ | $\sigma_{z}/\sigma_{R}$ |  $M_{\bullet}$ ($M_{\odot}$) |  $\Upsilon_{i^{\prime}}$ $(M_{\odot}/L_{\odot,i^{\prime}})$
20 | Spherical | Kept | $73^{+12}_{-14}$ | $1.0^{+0.3}_{-0.2}$ | - | $(2.8^{+0.8}_{-0.8})\times 10^{10}$ | $4.9^{+0.9}_{-0.7}$
20 | Spherical | Removed | $71^{+13}_{-15}$ | $0.83^{+0.23}_{-0.18}$ | - | $(2.0^{+0.9}_{-0.9})\times 10^{10}$ | $4.7^{+0.8}_{-0.7}$
20 | Cylindrical | Kept | $66^{+15}_{-13}$ | - | $0.92^{+0.09}_{-0.10}$ | $(2.7^{+0.5}_{-0.4})\times 10^{10}$ | $4.8^{+0.5}_{-0.5}$
20 | Cylindrical | Removed | $66^{+17}_{-15}$ | - | $0.91^{+0.08}_{-0.10}$ | $(2.4^{+0.5}_{-0.5})\times 10^{10}$ | $4.9^{+0.5}_{-0.5}$
18 | Spherical | Kept | $72^{+12}_{-14}$ | $0.85^{+0.20}_{-0.16}$ | - | $(2.0^{+0.6}_{-0.6})\times 10^{10}$ | $4.1^{+0.7}_{-0.5}$
18 | Spherical | Removed | $72^{+12}_{-15}$ | $0.83^{+0.18}_{-0.16}$ | - | $(1.9^{+0.6}_{-0.7})\times 10^{10}$ | $4.1^{+0.6}_{-0.5}$
18 | Cylindrical | Kept | $66^{+15}_{-15}$ | - | $0.89^{+0.08}_{-0.10}$ | $(2.3^{+0.3}_{-0.3})\times 10^{10}$ | $4.2^{+0.4}_{-0.4}$
18 | Cylindrical | Removed | $68^{+15}_{-16}$ | - | $0.89^{+0.07}_{-0.09}$ | $(2.2^{+0.4}_{-0.4})\times 10^{10}$ | $4.2^{+0.4}_{-0.4}$
Table 4: The Self-Consistent Model posterior medians and 1-$\sigma$ ranges for
the different kinematic models/data sets.

###  3.2 Results

We perform MCMC fitting of the JAM models, with a $sin(i)$ prior on the
inclination and wide uniform priors on the other parameters. The marginal
posterior values and 68% confidence regions are summarized in Table 4 for two
S/N binnings, both velocity ellipsoid models, and retention or exclusion of
the central kinematic point. We note that while error bars are larger on the
spherical model parameters and while low S/N binning tends to slightly reduce
the best-fit black hole masses, all results are consistent at the $\sim
1\sigma$ level.

Figure 4: $V_{\rm rms}=\sqrt{V^{2}+\sigma^{2}}$ $vs.$ radius compared with the
best-fit spherical JAM model. While the fit was to all data, the points are
color-coded by major axis/minor axis sector (Fig. 1), and the curves show the
expected model runs along the major and minor axes. The $V$ and $\sigma$
values are available as a DbF file.

Some additional systematic sensitivity comes from the choice of the
spectrum/template polynomials. We explored this by re-fitting the cylindrical
$\text{S/N}=20$ model (with its smaller statistical errors), varying the
degree of multiplicative and additive polynomials from 4 through 7 (16
combinations). Since the SMBHB mass $M_{\bullet}$ is our prime measured
parameter, we focus on its sensitivity. While over half of the values fell
within $0.5\times 10^{10}M_{\odot}$ of the 4A+6M fit result, values ranged
from $2.25-3.75\times 10^{10}M_{\odot}$. Evidently, significant systematic
uncertainties remain in the velocity extractions; these can be controlled with
higher S/N data.

Since 4C+37.11 appears to be a slow rotator, we take the spherical S/N=20 fits
as our standard, giving $M_{\bullet}=(2.8\pm 0.8)\times 10^{10}$ $M_{\odot}$,
statistical errors only. While the mass is large, so are the error bars,
encompassing nearly all of the other fits and polynomial choices. However, a
slightly more conservative interpretation might adopt $M_{\bullet}\approx
2.4\times 10^{10}M_{\odot}$, about $\sigma/2$ lower, since the central S/N=20
bin does tend to increase the mass. The radial $V_{\rm rms}$ run and corner
plots for the spherical fit are shown in Figs. 4 and 5. In the latter we note
that $i$ is only weakly constrained within the prior and that $M_{\bullet}$ is
correlated with the velocity anisotropy, with larger hole masses demanded for
$\beta<0$. Higher S/N is clearly needed to pin down $\beta$ (and its possible
radial dependence due to core scouring, see below).

Figure 5: Corner plot for the Spherical S/N=20 JAM fit.

Figure 6: (a) Correlation between core radius and black hole mass from Dullo
et al. (2021), including measurement of depleted core Holm 15A from Mehrgan et
al. (2019). (b) Correlation with effective velocity dispersion for BCGs from
Bogdán et al. (2018), represented by the orange data points and blue fit line.
The black line is the trend for all core ellipticals from Saglia et al.
(2016). Note that 4C+37.11/Holm 15 lie well above the relation for all
ellipticals, but are consistent with the BCG trend.

Using the black hole mass, constant stellar mass-to-light ratio, and
inclination for the spherical model, we estimate the radius of the sphere of
influence $r_{SOI}$ of the binary. Following Mehrgan et al. (2019) and Simon
et al. (2024), we define $r_{SOI}$ by $M_{\star}(\leq r_{SOI})=M_{\bullet}$,
where $M_{\star}$ is the enclosed stellar mass. Using the routine
mge_radial_mass of JamPy (Cappellari, 2008) to compute the enclosed stellar
mass within a radius $r$ using the MGE of the surface brightness, we find
$r_{SOI}\approx 1.29^{\prime\prime}$ ($1.38$ kpc), quite consistent with the
observed $V_{rms}$ flattening.

##  4 Discussion and Conclusions

It is interesting to compare our dynamical $M_{\bullet}\approx 2.8\times
10^{10}M_{\odot}$ estimate for 4C+37.11’s binary with values found for other
massive ellipticals. These $M_{\bullet}$ values correlate with host
properties, especially those of the host core, giving insight into black
holes/host co-evolution (Kormendy & Ho, 2013; Saglia et al., 2016; Bogdán et
al., 2018).

The core radius $r_{b}$ has been found to tightly correlate to the central
black hole mass among massive ellipticals (Mehrgan et al., 2019; Dullo et al.,
2021). We compare our 4C+37.11 values with dynamically measured $M_{\bullet}$
and $r_{b}$ for other elliptical galaxies (Dullo et al., 2021), including Holm
15A, which has the largest dynamically measured black hole mass in the local
universe, from Mehrgan et al. (2019). We find our value is quite consistent
with these authors’ best-fit trend for cored ellipticals (Figure 6a).

Another well-known correlation is the $M_{\bullet}-\sigma$ relation, which has
been found to depend on the host type. While our $\sigma_{e}=332\pm 11$ km/s
is somewhat suspect since we can only integrate to $\approx 0.25R_{e}$, we can
compare with the trends seen by various authors. We first consider the cored
ellipticals with dynamical mass measurements in Saglia et al. (2016). Our
estimate exceeds their best-fit correlation’s prediction by $\sim 7.2\times$.
However, if we compare with the cored ellipticals which are also brightest
cluster galaxies (BCGs) from Bogdán et al. (2018), we find that 4C+37.11 (and
Holm 15A) follow their best-fit trend more closely (Figure 6b).

The $M_{\bullet}\approx 2.8\times 10^{10}M_{\odot}$ dynamical mass for
4C+37.11’s SMBHB is one of the largest measured in the local universe. While
exceeding the value predicted for a typical elliptical, it is in line with
expectations for cored ellipticals and, especially, those that are also BCGs.
This fits well with the picture that the host is a fossil cluster, the product
of several major mergers. If these mergers were largely dry, dissipation-less
events, they would grow the central black hole mass faster than the stellar
velocity dispersion (Mehrgan et al., 2019). In addition, the back-action from
the current binary (and likely from earlier binary phases) has ‘scoured’ the
core, removing stars capable of exerting dynamical friction via three-body
gravitational slingshots (Begelman et al., 1980). In that, 4C+37.11 resembles
other extreme cored ellipticals, including Holm 15A. This picture can be
tested and extended by additional kinematic studies, which can measure
$\beta(r)$, tighten the $M_{\bullet}$ estimate and constrain the orbit
anisotropy of the stars in the host core, probing the system’s merger history.
While such measurements will be challenging or impossible from the ground, the
JWST NIRSpec IFU field is well matched to the host $r_{b}$ and can provide
important information on the core kinematics.

We thank the anonymous referee for constructive comments which substantially
improved the paper and thank the Stanford on the Moon program for publication
support. AP acknowledges the National Science Foundation for providing support
to staff for independent research. This paper used data from program
GN-2015B-FT-6, which was recovered from the Gemini Observatory Archive and
processed with the Gemini IRAF package. Thus this work is based in part on
observations obtained at the international Gemini Observatory, a program of
NSF’s NOIRLab, which is managed by the Association of Universities for
Research in Astronomy (AURA) under a cooperative agreement with the National
Science Foundation on behalf of the Gemini Observatory partnership: the
National Science Foundation (United States), National Research Council
(Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio
de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência,
Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space
Science Institute (Republic of Korea).

## References

  * Andrade-Santos et al. (2016) Andrade-Santos, F., Bogdán, Á., Romani, R. W., et al. 2016, ApJ, 826, 91, doi: 10.3847/0004-637X/826/1/91
  * Bansal et al. (2017) Bansal, K., Taylor, G. B., Peck, A. B., Zavala, R. T., & Romani, R. W. 2017, ApJ, 843, 14, doi: 10.3847/1538-4357/aa74e1
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307, doi: 10.1038/287307a0
  * Bennett et al. (2014) Bennett, C. L., Larson, D., Weiland, J. L., & Hinshaw, G. 2014, ApJ, 794, 135, doi: 10.1088/0004-637X/794/2/135
  * Bogdán et al. (2018) Bogdán, Á., Lovisari, L., Volonteri, M., & Dubois, Y. 2018, ApJ, 852, 131, doi: 10.3847/1538-4357/aa9ab5
  * Bradley et al. (2023) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2023, astropy/photutils: 1.8.0, 1.8.0, Zenodo, doi: 10.5281/zenodo.7946442
  * Cappellari (2002) Cappellari, M. 2002, MNRAS, 333, 400, doi: 10.1046/j.1365-8711.2002.05412.x
  * Cappellari (2008) —. 2008, MNRAS, 390, 71, doi: 10.1111/j.1365-2966.2008.13754.x
  * Cappellari (2016) Cappellari, M. 2016, Annual Review of Astronomy and Astrophysics, 54, 597, doi: 10.1146/annurev-astro-082214-122432
  * Cappellari (2017) —. 2017, Monthly Notices of the Royal Astronomical Society, 466, 798, doi: 10.1093/mnras/stw3020
  * Cappellari (2023) Cappellari, M. 2023, MNRAS, 526, 3273, doi: 10.1093/mnras/stad2597
  * Cappellari & Copin (2003) Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345, doi: 10.1046/j.1365-8711.2003.06541.x
  * Cappellari et al. (2010) Cappellari, M., McDermid, R. M., Bacon, R., et al. 2010, in AIP Conference Proceedings (AIP), doi: 10.1063/1.3458489
  * Cappellari et al. (2013) Cappellari, M., Scott, N., Alatalo, K., et al. 2013, Monthly Notices of the Royal Astronomical Society, 432, 1709, doi: 10.1093/mnras/stt562
  * Dullo et al. (2021) Dullo, B. T., Gil de Paz, A., & Knapen, J. H. 2021, ApJ, 908, 134, doi: 10.3847/1538-4357/abceae
  * Erwin (2015) Erwin, P. 2015, The Astrophysical Journal, 799, 226, doi: 10.1088/0004-637x/799/2/226
  * Graham et al. (2003) Graham, A. W., Erwin, P., Trujillo, I., & Ramos, A. A. 2003, The Astronomical Journal, 125, 2951, doi: 10.1086/375320
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, Annual Review of Astronomy and Astrophysics, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Maness et al. (2004) Maness, H. L., Taylor, G. B., Zavala, R. T., Peck, A. B., & Pollack, L. K. 2004, The Astrophysical Journal, 602, 123, doi: 10.1086/380919
  * Mehrgan et al. (2019) Mehrgan, K., Thomas, J., Saglia, R., et al. 2019, The Astrophysical Journal, 887, 195, doi: 10.3847/1538-4357/ab5856
  * Rodriguez et al. (2006) Rodriguez, C., Taylor, G. B., Zavala, R. T., et al. 2006, ApJ, 646, 49, doi: 10.1086/504825
  * Romani et al. (2014) Romani, R. W., Forman, W. R., Jones, C., et al. 2014, ApJ, 780, 149, doi: 10.1088/0004-637X/780/2/149
  * Saglia et al. (2016) Saglia, R. P., Opitsch, M., Erwin, P., et al. 2016, ApJ, 818, 47, doi: 10.3847/0004-637X/818/1/47
  * Simon et al. (2024) Simon, D. A., Cappellari, M., & Hartke, J. 2024, MNRAS, 527, 2341, doi: 10.1093/mnras/stad3309
  * Trujillo et al. (2004) Trujillo, I., Erwin, P., Asensio Ramos, A., & Graham, A. W. 2004, AJ, 127, 1917, doi: 10.1086/382712
  * Valdes et al. (2004) Valdes, F., Gupta, R., Rose, J. A., Singh, H. P., & Bell, D. J. 2004, The Astrophysical Journal Supplement Series, 152, 251, doi: 10.1086/386343


# Generalist Multimodal AI: A Review of Architectures, Challenges and
Opportunities

Sai Munikoti Pacific Northwest National LabRichlandWAUSA<EMAIL_ADDRESS>,  Ian Stewart Pacific Northwest National LabRichlandWAUSA
<EMAIL_ADDRESS>,  Sameera Horawalavithana Pacific Northwest National
LabRichlandWAUSA<EMAIL_ADDRESS>,  Henry Kvinge Pacific
Northwest National LabRichlandWAUSA<EMAIL_ADDRESS>,  Tegan Emerson
Pacific Northwest National LabRichlandWAUSA<EMAIL_ADDRESS>and  Karl
Pazdernik Pacific Northwest National LabRichlandWAUSA<EMAIL_ADDRESS>
###### Abstract.

Multimodal models, which learn from two or more modalities, Multimodal models
are expected to be a critical component to future advances in artificial
intelligence. This field is starting to grow rapidly with a surge of new
design elements motivated by the success of foundation models in NLP natural
language processing (NLP) and vision. It is widely hoped that further
extending the foundation models to multiple other modalities modalities (e.g.,
text, image, video, sensor, time series, graph, etc.) will ultimately lead to
artificial generalist multimodal models, i.e. one model across different data
modalities and tasks. However, there is little research that systematically
analyzes recent trends in the multimodal literaturemultimodal models
(particularly the ones that work beyond text and vision) with respect to the
underling architecture proposed. Therefore, this work provides a fresh
perspective on foundation generalist multimodal models (FMMsGMMs) via a novel
architecture and training configuration specific taxonomy. This includes
factors such as Unifiability, Modularity, and Adaptability that are pertinent
and essential to the wide adoption and application of FMMsGMMs. The review
further highlights key challenges and prospects for the field and guide the
researchers into the new advancements.

multimodal, transformer, multitask, unified architectures, zero shot


# Quantum optimization of coherent chaotic systems: A case for buses of
Kathmandu.

Kiran Adhikari  Institute for Communications Engineering, Technical University
of Munich, Arcisstraße 21, 80333 München  Aman Ganeju  Khwopa College,
Tribhuvan University, Kathmandu, Nepal  Iva Kumari Lamichhane  Tri-Chandra
Multiple Campus, Tribhuvan University, Kathmandu, Nepal  Rohit Bhattarai  Tri-
Chandra Multiple Campus, Tribhuvan University, Kathmandu, Nepal  Manghang
Limbu  Tri-Chandra Multiple Campus, Tribhuvan University, Kathmandu, Nepal
Nishma Bhattarai  St. Xavier’s College, Maitighar, Kathmandu, Nepal  Christian
Deppe  Institute for Communications Engineering, Technical University of
Munich, Arcisstraße 21, 80333 München

###### Abstract

In this paper, we propose a novel quantum computing approach to solve the
real-world problem of optimizing transportation in bustling Kathmandu city.
The transportation system in Kathmandu is chaotic, with no central authority
controlling the transportation. We leverage this chaotic feature in our
quantum optimization procedure. The quantum chaos theory’s Wigner-Dyson
distribution surfaced as the most effective bus spacing distribution for a bus
driver to maximize their profit. We investigate the statistical properties of
the buses with real-time GPS bus location data and optimize bus spacing and
interval distribution around the 27 km circular ring road in Kathmandu. Using
tools like quantum simulation, eigenvalue distributions, and output wave
function analysis, we show that such optimal bus spacing distribution could be
achieved.

##  1 Introduction

In the early 1980s, Benioff and Feynman introduced the concept of quantum
computing, highlighting that quantum mechanical computers outperformed
classical ones for specific tasks [1]. Feynman suggested applying quantum
mechanics to computational problems. As stated in [2], it was suggested that
standard quantum systems could simulate complex quantum systems to solve
problems beyond the capabilities of classical computers. Although the
realization of such a quantum simulator was unclear at the time, Feynman’s
concept significantly influenced the development of quantum computing.
Subsequently, the concept of a quantum Turing machine was introduced, and the
existence of universal models based on quantum mechanics was theoretically
established [3]. In simple terms, quantum models can achieve the same
computational capabilities as classical computers. The first quantum algorithm
was presented in [4], which demonstrated exponential acceleration in quantum
computation compared to classical computation. Following that, quantum
algorithms were proposed in [5] and [6], highlighting the advantages of
quantum computing in solving specific problems. However, these problems were
rather artificial and their practical applicability is limited.

In 1994, Shor introduced a quantum algorithm (named after him today), designed
for factorizing large integers [7]. Large prime factorization is known to be
an NP-hard problem and is the foundation for the security of the RSA public
key cryptosystem. Solving this problem using classical computers requires
exponential time. However, the Shor algorithm demonstrated that quantum
computing could achieve the same task in polynomial time, posing a potential
threat to the security of RSA encryption by making it easier to decipher.
Subsequently, Grover proposed a quantum searching algorithm that could
efficiently locate a specific item in an unordered database. The Grover
algorithm takes advantage of quantum parallelism, enabling it to check all the
data simultaneously in each iteration. This feature significantly reduces the
computational complexity of solving the search problem. The unique computation
capabilities of Shor and Grover algorithms have since captured widespread
attention and interest. One of the areas where quantum algorithms provide
substantial improvements is the problem of optimization [8, 9, 10, 11, 12, 13,
14, 15, 16].

Typically, optimization problems involve finding the extremum (maximum or
minimum) of an objective function, which could represent various quantities
like cost, energy, profit, or time. In many cases, this function depends on a
substantial number of parameters, making a brute-force approach impractical
due to the need to search through a high-dimensional parameter space. Quantum
computers operate in Hilbert space, and this space grows exponentially with
the number of qubits ($2^{N}$). The concept behind quantum computing is to
leverage this vast state space through quantum entanglement to enhance the
efficiency of finding the optimal solution, ideally achieving exponential
speed-up [17, 18, 19, 20, 10]. However, upon closer examination, it becomes
apparent that the speed-up for classical optimization problems is often quite
modest[21, 22, 23]. In contrast, quantum speed-up becomes most advantageous
for problems that are directly tied to the quantum-mechanical description of
natural phenomena. An excellent example is finding the many-electron wave
function of a molecular system. Classical computers struggle to solve such
problems accurately when the number of electrons exceeds a few tens due to the
exponential growth of Hilbert space with electron count. Quantum computers,
with their large state space, can simulate chemical systems and compute their
properties, including correlations and reaction rates. Of course, this is
under the condition that the challenge of efficiently mapping the fermionic
problem to the available qubit hardware is overcome.

Most of our currently best-known quantum algorithms exploit huge structures in
the problem to solve it exponentially faster than classical computers [24].
However, the history of physics has taught us that tools from physics are
remarkably powerful even when there is little structure in the underlying
system. A compelling exemplification of this principle is found in the domains
of thermodynamics and statistical mechanics. In this paper, we look for
quantum algorithms for a chaotic system where there is very little structure
in the underlying system.

Random matrix theory presents itself as a promising domain to look for such a
problem. Random-matrix theory deals with the statistical properties of large
matrices with randomly distributed elements. The probability distribution of
the matrices is taken as input, from which the correlation functions of
eigenvalues and eigenvectors are derived as output. From the correlation
functions one then computes the physical properties of the system. In the
context of systems being well modeled by the Random Matrix Theory, the bus
system in Cuernavaca, Mexico [25] in the late 1990s has become a canonical
physical system [26, 27, 28, 29]. Moreover, since the publishing of the
seminal paper, the associated phenomenon is still being rediscovered in
different systems like cars on a motorway [30] and subway trains in New York
City [31]. A similar phenomenon can be expected in the transportation system
of Kathmandu. It is with this motivation we look for quantum optimization of
transportation systems.

With rapid urbanization, transportation supply and demand imbalances have led
to common problems like pollution, traffic congestion, and energy waste in
various countries. The best solution is prioritizing public transport and
improving urban public transport systems. But the problem is that the urban
transit network design problem (UTNDP) which is concerned with the
determination of a set of routes with corresponding schedules for an urban
public transport system is an NP-Hard problem which means that in realistic
time frames, it is not possible to find an optimal verifiable solution. [32]
The large number of possible solutions and computational expenses make it
impractical to explore the entire solution space. Therefore, a lot of
different methods have been devised to tackle the UTNDP. (For a more detailed
overview of these methods, see [32] and references therein).

The UTNDP can be subdivided into two major components, namely the transit
routing problem and the transit scheduling problem [33]. Generally, the urban
transit routing problem involves the development of efficient transit routes
(e.g., bus routes) on an existing road network, with predefined pickup/drop-
off points (e.g., bus stops). On the other hand, the urban transit scheduling
problem (UTSP) is charged with assigning the schedules for the passenger-
carrying vehicles. Traditionally, the two phases are usually implemented
sequentially (or iteratively), with the routes determined in advance of the
schedule. In this paper, we concentrate on the urban transit routing problem
and present a novel approach by giving an algorithm to model the problem in a
quantum mechanical framework to make it more amenable to further quantum
processing.

The structure of this paper is as follows. In section 2, we discuss the basics
of Random Matrix Theory, and in section 2.1, we discuss the transport system
of Cuernavaca, which Random Matrix Theory can model. In Section 3, we model
the bus system of Kathmandu as a coherent chaotic system using the Dyson gas
model. In section 4, we suggest an algorithm for optimizing the route and also
suggest its relation to signatures of quantum chaos. In section 5, we give an
outlook for future directions, while in section 6, we summarize our paper.

##  2 Random matrix theory

It is believed that the dynamics of general many-body quantum chaotic systems
beyond a timescale $t_{Th}$, the so-called “Thouless time,” follow predictions
from random matrix theory (RMT) [34] i.e., their late-time behavior resembles
that of a random matrix chosen from an ensemble consistent with the system’s
symmetries [35, 36, 37].

Random matrix theory is used to characterize complex quantum systems when
there is limited knowledge about the Hamiltonian. The fundamental hypothesis
is that the Hamiltonian can be treated as a random matrix drawn from an
ensemble with appropriate symmetries [38]. This approach is particularly
useful for systems with many degrees of freedom and unknown interaction
couplings among them.

Within this approach, one is often interested in the eigenspectrum of the
Hamiltonian. Let the eigenvalues be $\lambda_{1}\leq\lambda_{2}\leq....$ and
let $S_{n}=\lambda_{n+1}-\lambda_{n}$ be the consecutive splittings. The
average value $D=\langle S_{n}\rangle$ is the mean splitting over the
eigenvalue intervals. In integrable Hamiltonian models, the eigenvalues are
uncorrelated, they are not prohibited from crossing and are usually described
by Poisson statistics [39]. In this case, the distribution $P(s)$ of the
neighboring spacing $s=S/D$, where $S$ is a particular spacing and $D$ is the
mean distance between neighboring intervals is given by:

$P(s)=\frac{1}{D}e^{-s}$ (1)

In contrast, the level spacing distribution P(s) of chaotic models is closely
approximated by the Wigner-Dyson (WD) distribution:

$P(s)=b_{\beta}s^{\beta}e^{-\alpha_{\beta}}s^{2}$ (2)

where $\beta$ depends on which universality class of random matrices the
chaotic Hamiltonian belongs to [40]: $\beta=1$ for the Gaussian Orthogonal
Ensemble (GOE), $\beta=2$ for the Gaussian Unitary Ensemble (GUE), and
$\beta=4$ for the Gaussian symplectic ensemble (GSE), and
$b_{1}=\pi/2,a_{1}=\pi/4$; $b_{2}=32/\pi^{2},a_{2}=4/\pi$;
$b_{4}=262144/729\pi^{3},a_{4}=64/9\pi$. To compare the spectral statistics
with regular and chaotic limits and also exhibit interpolation between them,
different distribution functions have been proposed. For example, for the
Gaussian Orthogonal Ensemble (GOE) statistics, one popular intermediate
distribution is the Broody distribution:

$\displaystyle
P(s)=b(1+q)s^{q}e^{-bs^{q+1}},b=\left(\Gamma\frac{2+q}{1+q}\right)^{q+1}$

where $q=0$ corresponds to the Poisson limit while $q=1$ corresponds to
Wigner-Dyson limit. Such distributions are shown in figure 1.

Figure 1: level spacing distribution $P(s)$ for various systems

The Bohigas-Giannoni-Schmit (BGS) conjecture [37], [41] associates the quantum
chaotic properties of a system with the correlations between its energy
levels. Chaotic Hamiltonian exhibits level correlations in agreement with the
predictions of random matrix theory (RMT) [40]. Adjacent eigenvalues show
level repulsion and, at larger energy scales, signals of spectral rigidity.
This understanding can be extended to the eigenspectrum of the Hamiltonian,
allowing us to determine the chaotic nature of a quantum state based on its
level spacing distribution [42].

We will study one of the public bus systems of Kathmandu, namely, the Ring
Road Bus system, and explore the chaotic nature of the system using various
quantum chaos diagnostics as a function of eigenvalue level spacing and wave-
function analysis.

###  2.1 Cuernavaca bus system and Dyson gas

Dyson gas model is a mathematical model proposed by Freeman Dyson to describe
the distribution of eigenvalues of random matrices in which Dyson gas is
mirrored as an interacting many-particle system moving in one dimension[41].

Buses of Cuernavaca are analogous to Dyson gas particles, thus the Dyson gas
model can describe the statistical properties like bus spacing distribution of
transport in Cuernavaca [25]. Since very few changes in the initial conditions
can bring out unpredictable changes in the Cuernavaca bus system as well as
Dyson gas behaviors, both behave like chaotic quantum systems. This kind of
complex system is characterized by minimum information content about the
Hamiltonian, as a result of which their statistical properties follow random
matrix ensembles( Gaussian Unitary Ensemble) prediction.

The motion of buses, as well as Dyson gas particles, is along one dimension
where $x_{1},x_{2},x_{3},...x_{n}$ denotes the position on the circle. Both
systems are interactive, where the interaction between buses is indirectly
facilitated by people assigned at different bus stops, which provide
information about distance and time gap with the preceding bus ($x_{1}-x_{2})$
and following bus ($x_{2}-x_{3}$). Likewise, the interaction potential that
exists between buses in Cuernavaca is the function of the spacing between the
buses [25]:

$V(x)\approx 1/x^{a}$ (3)

In a similar sense, the interaction potential between Dyson gas particles is
also a function of the distance between them which can be expressed as:

$V(x)\approx 1/\lambda^{a}$ (4)

where V(x) is interaction potential.

##  3 Modeling buses of Kathmandu as a coherent chaotic system

It is well known that a probability distribution of the distances between
subsequent buses with a centralized bus system is close to the Poisson
distribution and can be described by the standard bus route model [43]. The
origin of Poisson distribution in the centralized bus system comes from the
fact that two neighboring buses are not interacting with each other. In
contrast to the centralized bus system like in other cities, there doesn’t
exist a centralized bus company in Kathmandu. Drivers work independently and
try to optimize the amount of money earned while following a given route.
Consequently, constraints such as a timetable that represents external
influence on transport do not exist. This leads to competition among the
drivers and to their mutual interaction which gives a departure from Poisson-
like distribution. A Poisson-like distribution implies, that the probability
of close encounters between two buses is high (bus clustering) which conflicts
with the effort of the drivers to maximize the number of transported
passengers and accordingly to maximize the distance to the preceding bus. This
on the other hand means, that they tend not to arrive at a given bus stop
immediately after their predecessors. To achieve that, they learn information
about the time of departure of the previous bus from collaborators stationed
at the bus stop, they arrive at. Knowing the time interval to the preceding
bus, the driver tries to optimize the distance to it by either slowing down or
speeding up. In such a way the obtained information leads to a direct
interaction between buses and changes the statistical properties.

The GUE properties of the bus arrival statistics can be well understood when
regarding the buses as a one-dimensional interacting gas. Buses here resemble
Dyson gas particles behavior and buses are one-dimensional gas particles
interacting with inverse square potential i.e. if $x$ is the distance between
the particles then $1/x^{2}$ is their interaction potential. It has already
been mentioned that the exact GUE statistics are obtained for a Coulomb
interaction between the gas particles, i.e., for the interaction potential $V$
given by:

$V=-\sum_{i<j}\log(|x_{i}-x_{j}|)+\frac{1}{2}\sum_{i}x_{i}^{2}$ (5)

In (5), the second term represents a force confining the gas close to the
origin. It is possible to consider more complicated Hamiltonian,

$H=\sum_{i=1}^{N}\frac{m_{i}v_{i}^{2}}{2}+\sum_{i=1}^{N}U(|x_{i-1}-x_{i}-l_{i}|)$
(6)

with $x_{i}$, $l_{i}$, $m_{i}$, and $v_{i}$ being the location, length, mass,
and velocity of the ith vehicle.

Other factors of the Kathmandu transportation system can also be considered
for realistic optimization. In some busy bus stops, there are traffic
policemen with some probability to check the bus’s stay time so that they can
minimize traffic jams. This can act as a cutoff for the stay time at those bus
stops. Drivers also communicate with buses that are returning in opposite
lanes to get information about other buses’ positions. This can be modeled by
considering a probabilistic interaction among buses that are not neighbors.
Every driver usually encounters the widespread phenomenon of so-called
bottleneck jam leading to transient clustering. So, one bus leaving earlier
may be caught by another bus behind. Furthermore, at a particular length
interval, there might be buses from different routes. In addition to the
inter-route competition, in this case, there will be competition among buses
from different routes too which we call, intra-route interaction. All of these
extra features can be considered for realistic optimization.

###  3.1 Ring Road route

In this paper, we will only consider one particular route of Kathmandu called
the Ring Road route. The circumference of the Ring Road route is 27 km with a
width of 31m on either side from the center of the road and was constructed in
1977 with China’s assistance [44].

The ring road of Kathmandu, as its name suggests, is a complete ring-like
transport system where the motion of vehicles is periodic along a closed path
around the cities of Kathmandu and Lalitpur. This closed circular path gives a
periodic boundary condition which makes it easier to do numerical and quantum
simulations. Furthermore, we can now think of buses as a one-dimensional gas
moving on a circle instead of a line, and we can safely ignore the second term
of Eq (5).

Figure 2: Map of Ring Road: The circumference of the Ring Road route is 27 km
with a width of 31m on either side from the center of the road.

Mahanagar Yatayat Samiti, a metropolitan transport committee, operates
approximately 55 buses along the Ring Road route. However, there are numerous
additional buses from different committees on the same route. It should be
remarked that these bus drivers still operate independently without adhering
to a fixed timetable, and thus cannot be modeled by standard bus route models
following Poisson distribution. A significant development by the committee is
the Mahanagar app, a smartphone application that enables users to track the
GPS positions of each bus. In figure 3, we have included time snapshots taken
from the Mahangar app. These snapshots contain the location of buses run by
the committee at various time stamps. In certain timestamps, such as in figure
4, the spacing between buses can be seen to be uniform as well. These kinds of
information, taken from the app, can be leveraged for optimization purposes.

(a) Clustering At 7:54 A.M

(b) Clustering At 8:30 A.M

(c) Clustering At 9:47 A.M

(d) Clustering At 10:30 A.M

Figure 3: Real-time GPS positions of approximately 55 buses from the Mahanagar
Yatayat Samiti, a metropolitan transport committee, as taken from the
Mahanagar app.  Figure 4: GPS position of buses at 9:21 AM, which shows an
example of approximately uniform spacing among buses.

###  3.2 What is to optimize?

The optimal spacing between buses can vary depending on whether you consider
the perspective of bus drivers or passengers. From the perspective of a
passenger waiting at a bus stop, the optimal spacing between buses would be
when the spacing has a probability distribution close to the Poisson
distribution. This type of distribution indicates a higher probability of
buses arriving at closer intervals, increasing the likelihood of frequent
encounters between buses. While this is better from the passenger’s
perspective, it is in conflict with the effort of the drivers to maximize the
number of transported passengers and, accordingly to maximize the distance to
the preceding bus. Therefore, from the driver’s perspective, the optimal
probability distribution for spacing between buses would be one that closely
approximates the Wigner-Dyson distribution. The third possibility would be a
balance between these two cases in which the probability distribution would be
close to the Broody distribution. Even though we have not provided proof that
these are indeed optimal scenarios, for now, we leave it as a conjecture.

In this paper, we will mostly focus on the bus driver’s perspective, who is
trying to maximize his profit while following a particular route. Somewhat
counter-intuitively, the optimal bus spacing distribution, in this case, would
be when the system is chaotic. In the next section, we will provide techniques
for such optimization procedures. The concept behind this approach involves
utilizing the real-time position data obtained from GPS systems installed in
buses. By employing a combination of quantum and classical algorithms and
simulations, we then generate recommendations for the optimal stopping time at
upcoming or current bus stops. These algorithms leverage the power of quantum
computing alongside traditional computing techniques to process and analyze
the bus position data, ultimately suggesting the most efficient and effective
stopping times for bus operations.

##  4 Optimizing the Ring road route

In this section, we will provide various quantum techniques to optimize the
ring road route. To simplify the model, we assume that buses travel in a
circular path, which is topologically equivalent to the ring road route.
Therefore, the principles of physics and optimization can be applied
similarly.

Let us denote the positions of N buses moving in a circle as
$x_{1},x_{2},...,x_{N}$, with the periodic boundary condition that
$x_{N+1}=x_{1}$. Let $G(x_{1},x_{2},...,x_{N})$ be the level spacing
distribution for these positions. In the ring road route, these real-time
positions of buses are given by GPS installed in buses. Now, the optimization
is to provide new positions $x_{i}$ such that the new level spacing
distribution is close to that of the Wigner-Dyson distribution for GUE which
we denote by $P_{\text{WD}}$. More concretely, if $x_{1},x_{2},...,x_{N}$ is
the position of the buses at time $t$, the task is to provide new bus
positions $x_{1}+\delta x_{1},x_{2}+\delta x_{2},....,x_{N}+\delta x_{N}$
which each bus should shift in time $t+\delta t$ such that

$|G(x_{1}+\delta x_{1},x_{2}+\delta x_{2},....,x_{N}+\delta
x_{N})-P_{\text{WD}}|\leq\epsilon$

.

There are two possible ways that one can shift the buses to the distribution
$G(x_{i})$ to $P_{WD}$:

  * •

Decrease or increase the speed of the buses depending on the sign of $\delta
x_{i}$

  * •

Manipulate stopping time at the next stop depending on $\delta x$

However, in a busy road system with speed limits on highways, it might not be
possible to use increase the speed of the buses. This leaves the option to
manipulate the stopping time of the next bus stop. To work out the algorithm
for it, we need to introduce a few relevant terms:

  1. 1.

Due to the non-uniform speed of the buses and the chaotic nature of traffic,
we would consider an average velocity between two stops $i$ and $j$ as
$\braket{v_{ij}}$ where $i<j$.

  2. 2.

Since a bus stops at a bus stop for a considerable amount of time every time,
we introduce a maximum allowed stopping time at each stop as $t_{m}$ and an
average stopping time at bus stops as $\braket{t_{s}}$. The maximum stopping
time $t_{m}$ may be due to city regulations that impose such limitations, or
it could be motivated by the practical consideration that excessively long
stopping times would allow buses from behind to catch up.

Depending on $\delta x$, a bus might need to stop at a bus stop for a longer
period of time or for a shorter period of time. To consider all the scenarios,
we considered three possible states:

#### Case I: $\delta x_{i}=0$

If $\delta x_{i}=0$ for a bus indexed $i$, then the level spacing distribution
resembles the WD distribution. The bus can stop at the next bus stop for an
average time $\braket{t_{s}}$.

#### Case II: $\delta x_{i}\neq 0$

A positive value of $\delta x_{i}$ means that the bus is lagging behind
considering the position given by the new optimized distribution and the
negative value means the bus is faster than the given position. Consequently,
the bus must stop at the bus stop for less or more time than $\braket{t_{s}}$
respectively. The bus would need to stop at the next bus stop for
$t=\braket{t_{s}}\pm\frac{\delta x_{i}}{\braket{v_{ij}}}$, where j is the
position of the next bus stop.

However, when $\delta x_{i}<0$, a bus cannot stop more than the maximum
stopping time $t_{m}$. If $t\geq t_{m}$, the bus must stop at the bus stop for
time $t_{m}$ and the remaining distance needed to lag,
$d=\braket{v_{ij}}(\braket{t_{s}}-t_{m})$ will be reconciled at the next bus
stop.

###  4.1 Spectral analysis

We start with the Hamiltonian,

$H=\sum_{i=1}^{N}\frac{m_{i}v_{i}^{2}}{2}+\sum_{i=1}^{N}U(|x_{i-1}-x_{i}|)$
(7)

with $x_{i}$, $m_{i}$, and $v_{i}$ being the location, mass, and velocity of
the ith vehicle with some guess for the potential function
$U(|x_{i-1}-x_{i}|)$. It could just be $|x_{i-1}-x_{i}|$ as a starting case.
Changing the parameters $x_{i}$ by $\delta x_{i}$ gives a new Hamiltonian
$H^{\prime}$. One can then check the eigen spectrum of this Hamiltonian using
quantum computer [45, 46]. If the level spacing distribution is close to the
Wigner-Dyson distribution, we will accept this new position of the buses.
Otherwise, we will keep changing the $\delta x_{i}$. We do this until we reach
the level spacing distribution close to the Wigner-Dyson distribution. When
this happens, the Hamiltonian is expected to be close to the Dyson gas
Hamiltonian with potential,

$U=-\sum_{i<j}\log(|x_{i}-x_{j}|)$ (8)

By solving this, one can obtain the new optimal parameters for bus positions.
In addition to the level spacing distribution analysis, we can use additional
diagnostics, using quantum computer , to check if chaos has emerged. Given the
level spacings $S_{i}=\lambda_{i+1}-\lambda_{i}$, we can define the following
ratios:

$r_{i}=\frac{\text{Min}(S_{i},S_{i+1})}{\text{Max}(S_{i},S_{i+1})}$ (9)

.

For Poisson-like level spacing distribution, $r_{i}\approx 0.38629$, while for
Wigner-Dyson-like level spacing distribution, $r_{i}\approx 0.60266$. The
obtained $r_{i}$ should be close to $0.60266$ to get the optimal bus spacing.
In this scheme, the task is reduced to obtaining the eigenvalues of the
Hamiltonian using quantum computer, analyzing it in classical machines, and
updating the Hamiltonian.

One approach to achieve this goal involves conducting a quantum simulation of
the Hamiltonian. By manipulating the parameters within the Hamiltonian, one
can analyze its eigenvalues spectrum. Alternatively, this simulation can be
directly performed on a quantum computer [47]. In this method, a unitary
$U=e^{-iHt}$ is simulated on a quantum computer. Utilizing quantum phase
estimation, eigenvalues can be obtained with a complexity of $O(\log N)^{2}$,
a significant improvement compared to the classical approach which requires
$O(N^{3})$ operations [47]. Similar to the variational technique, one have the
flexibility to adjust both the parameters within the Hamiltonian and the
unitary operator, allowing for a more adaptable and efficient execution of the
algorithm.

###  4.2 Output wave function analysis

In the previous section, we considered analyzing the output spectrum for the
optimization process. However, the technique described therein involved
complicated quantum circuits, like quantum phase estimation with significant
depth, making it difficult to implement in the near term. One could also
directly analyze the output wave function generated by a parameterized quantum
circuit, which can be performed efficiently within polynomial time [48, 49].

One potential metric for analysis is the survival probability (SP), defined as
the overlap between a given initial state $\ket{\psi}$ and the output quantum
state $\ket{\phi}$:

$\text{SP}=|\langle\psi|\phi\rangle|^{2}$

This quantity indicates the likelihood of finding the initial system in its
original state even after undergoing a unitary transformation. Employing the
Hadamard test, one can calculate SP. A value of $1$ signifies that the output
state closely matches the initial state, while $0$ indicates complete
orthogonality. In the context of a chaotic quantum system, one would expect SP
to be close to $0$. Consequently, the Hamiltonian is adjusted iteratively
until SP decreases, and the corresponding Hamiltonian parameters are noted.

It is crucial to note that these parameters might not necessarily be optimal.
Achieving true optimality would necessitate a comprehensive spectral analysis,
a task demanding powerful quantum computers. Nevertheless, valuable insights
can still be extracted from the partial information gathered, guiding the
optimization process. Additionally, alternative metrics such as inverse
participation ratio, spread complexity, and Out of Time-Ordered correlation
experiments [50] can also be employed for analyzing the output wave function.

##  5 Outlook

Our research aims to state the conditions in which we can claim that the
Kathmandu bus system is indeed chaotic via the exploration of quantum
signatures. To identify these signatures, it is essential to gather
comprehensive data on the system, including bus positions, passenger loads,
traffic conditions, schedules, and other relevant parameters. The study’s
exploration of bus systems as coherent chaotic systems opens the door to more
detailed traffic flow analysis. Researchers can delve into the dynamics of
traffic which display critical points, non-equilibrium phase transitions,
noise-induced transitions, and fluctuation-induced ordering phenomena. By
adopting a quantum many-body system perspective, researches can uncover
quantum-like behavior in complex traffic patterns and contribute to the field
of traffic science.

Another research direction is to validate the proposed quantum optimization
paradigms using real data. By applying classical and quantum optimization
algorithms to this real data, a direct comparison can be made to assess the
effectiveness of both approaches in improving bus system performance. Of
course, it is not feasible to actually run the algorithm on a quantum computer
due to the severe limitations of quantum computers today. But the algorithm we
suggest is to optimize the bus system to achieve Wigner-Dyson distribution
which is observed a lot in quantum systems. This quantum-inspired algorithm
can be tested against classical optimization techniques and it is interesting
to compare the differences between them.

Besides, future research can also take a holistic approach by considering
optimization from the perspective of transportation authorities as well who
seek to minimize congestion and environmental impact unlike passengers who may
value shorter waiting times and drivers who aim to maximize profit. It is
interesting to explore multi-objective optimization via quantum route and to
see if these perspectives can lead to inspirations for more robust and
adaptable bus system designs.

In summary, the future research in this area is promising and multi-faceted.
Simulating Dyson gas, conducting in-depth traffic flow analysis, collecting
real data for validation, and approaching optimization from various
perspectives are all avenues that can further advance our understanding of bus
systems and potentially lead to more efficient and sustainable urban
transportation solutions.

##  6 Conclusion

In an endeavor to optimize the transportation system along Kathmandu’s Ring
Road, this study has investigated a novel approach integrating principles from
quantum chaos and Random Matrix theory. Kathmandu’s Ring Road serves as a
perfect model for optimization due to its unique, closed-loop design allowing
for a simplified model of buses moving in a one-dimensional gas format,
enabling numerical and quantum simulations. The road sees approximately 55
buses from the Mahanagar Yatayat Samiti, among other vehicles, which often
operate without a fixed timetable. The study explored three optimization
scenarios: one from the passenger’s perspective, another from the bus driver’s
view, and the third being a balance between the two. From a passenger’s
standpoint, ideal bus intervals would follow a Poisson distribution, ensuring
frequent bus arrivals. For drivers, the Wigner-Dyson distribution was deemed
optimal as they aim to maximize the distance to the preceding bus. The paper
primarily leaned towards the bus driver’s perspective. A smartphone
application called the Mahanagar app, which tracks the GPS locations of the
buses in real-time, was used to study the properties of the bus system. The
suggested technique involves adjusting the positions of the buses so that
their spacing closely aligns with the desired Wigner-Dyson distribution. Two
practical methods to achieve the desired spacing were proposed: manipulating
the bus speeds or adjusting their stopping times at the bus stops. Given road
conditions and regulatory restrictions, the latter was deemed more feasible.
Quantum algorithms, combined with classical methods, were suggested to analyze
bus position data and generate recommendations for optimal stopping times. The
Hamiltonian was shown to be able to be formulated with variables representing
the positions of vehicles or their spacing. By adjusting the parameters of the
Hamiltonian, the study aimed to bring its eigen spectrum closer to the Wigner-
Dyson distribution. A secondary method was proposed for scenarios where deep
quantum circuits, like those in quantum phase estimation, are infeasible. This
involves analyzing the output wave function of a parameterized circuit and
computing metrics like the survival probability using the SWAP test. This
provides insights into how chaotic the quantum system is, which could help in
the optimization process. In conclusion, by suggesting application of quantum
principles with real-world traffic data, this paper has paved the way for
innovative solutions to optimize the flow of buses on Kathmandu’s Ring Road.
While full implementation requires powerful quantum computers, the insights
garnered from this approach can provide significant improvements to the
current transport system.

##  7 Acknowledgements

This research is part of the Munich Quantum Valley, which is supported by the
Bavarian state government with funds from the Hightech Agenda Bayern Plus. We
also acknowledge support from Physics Without Frontiers (PWF) program of the
International Centre for Theoretical Physics (ICTP), Italy.

## References

  * [1] Paul Benioff.  “The computer as a physical system: A microscopic quantum mechanical hamiltonian model of computers as represented by turing machines” (1980).
  * [2] Richard P Feynman et al.  “Simulating physics with computers” (2018).
  * [3] David Deutsch.  “Quantum theory, the church–turing principle and the universal quantum computer” (1985).
  * [4] David Deutsch and Richard Jozsa.  “Rapid solution of problems by quantum computation” (1992).
  * [5] Ethan Bernstein and Umesh Vazirani.  “Quantum complexity theory” (1993).
  * [6] Daniel R Simon.  “On the power of quantum computation” (1997).
  * [7] Peter W Shor.  “Algorithms for quantum computation: discrete logarithms and factoring” (1994).
  * [8] Christoph Dürr, Mark Heiligman, Peter Hoyer, and Mehdi Mhalla.  “Quantum query complexity of some graph problems” (2006).
  * [9] Stephen P Jordan.  “Fast quantum algorithm for numerical gradient estimation” (2005).
  * [10] Aram W Harrow, Avinatan Hassidim, and Seth Lloyd.  “Quantum algorithm for linear systems of equations” (2009).
  * [11] Andrew M Childs, Robin Kothari, and Rolando D Somma.  “Quantum algorithm for systems of linear equations with exponentially improved dependence on precision” (2017).
  * [12] Fernando GSL Brandao and Krysta M Svore.  “Quantum speed-ups for solving semidefinite programs” (2017).
  * [13] Joran Van Apeldoorn, András Gilyén, Sander Gribling, and Ronald de Wolf.  “Quantum sdp-solvers: Better upper and lower bounds” (2020).
  * [14] Iordanis Kerenidis and Anupam Prakash.  “A quantum interior point method for lps and sdps” (2020).
  * [15] Joran van Apeldoorn, András Gilyén, Sander Gribling, and Ronald de Wolf.  “Convex optimization using quantum oracles” (2020).
  * [16] Shouvanik Chakrabarti, Andrew M Childs, Tongyang Li, and Xiaodi Wu.  “Quantum algorithms and lower bounds for convex optimization” (2020).
  * [17] Seth Lloyd.  “Universal quantum simulators” (1996).
  * [18] Daniel S Abrams and Seth Lloyd.  “Simulation of many-body fermi systems on a universal quantum computer” (1997).
  * [19] Daniel S Abrams and Seth Lloyd.  “Quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors” (1999).
  * [20] Alán Aspuru-Guzik, Anthony D Dutoi, Peter J Love, and Martin Head-Gordon.  “Simulated quantum computation of molecular energies” (2005).
  * [21] Vasil S Denchev, Sergio Boixo, Sergei V Isakov, Nan Ding, Ryan Babbush, Vadim Smelyanskiy, John Martinis, and Hartmut Neven.  “What is the computational value of finite-range tunneling?” (2016).
  * [22] Tameem Albash and Daniel A Lidar.  “Demonstration of a scaling advantage for a quantum annealer over simulated annealing” (2018).
  * [23] John A Smolin and Graeme Smith.  “Classical signature of quantum annealing” (2014).
  * [24] Scott Aaronson.  “How much structure is needed for huge quantum speedups?” (2022).
  * [25] Milan Krbálek and Petr Seba.  “The statistical properties of the city transport in cuernavaca (mexico) and random matrix ensembles” (2000).
  * [26] Milan Krbalek, Petr Šeba, and Peter Wagner.  “Headways in traffic flow: Remarks from a physical perspective” (2001).
  * [27] Milan Krbalek and Petr Seba.  “Headway statistics of public transport in mexican cities” (2002).
  * [28] Milan Krbálek.  “Inter-vehicle gap statistics on signal-controlled crossroads” (2008).
  * [29] Gernot Akemann, Jinho Baik, and Philippe Di Francesco.  “The oxford handbook of random matrix theory” (2011).
  * [30] AY Abul-Magd.  “Modeling highway-traffic headway distributions using superstatistics” (2007).
  * [31] Aukosh Jagannath and Thomas Trogdon.  “Random matrices and the new york city subway system” (2017).
  * [32] Lang Fan and Christine L Mumford.  “A metaheuristic approach to the urban transit routing problem” (2010).
  * [33] Partha Chakroborty.  “Genetic algorithms for optimal urban transit network design” (2003).
  * [34] D. J. Thouless.  “Maximum metallic resistance in thin wires” (1977).
  * [35] Jordan Cotler, Nicholas Hunter-Jones, Junyu Liu, and Beni Yoshida.  “Chaos, complexity, and random matrices” (2017).
  * [36] Freeman J Dyson.  “A brownian-motion model for the eigenvalues of a random matrix” (1962).
  * [37] Oriol Bohigas, Marie-Joya Giannoni, and Charles Schmit.  “Characterization of chaotic quantum spectra and universality of level fluctuation laws” (1984).
  * [38] Eugene P Wigner.  “On the statistical distribution of the widths and spacings of nuclear resonance levels” (1951).
  * [39] Michael Victor Berry and Michael Tabor.  “Level clustering in the regular spectrum” (1977).
  * [40] Madan Lal Mehta.  “Random matrices” (2004).
  * [41] Thomas Guhr, Axel Müller-Groeling, and Hans A Weidenmüller.  “Random-matrix theories in quantum physics: common concepts” (1998).
  * [42] Xiao Chen and Andreas WW Ludwig.  “Universal spectral correlations in the chaotic wave function and the development of quantum chaos” (2018).
  * [43] OJ O’loan, Martin R Evans, and Michael E Cates.  “Jamming transition in a homogeneous one-dimensional system: The bus route model” (1998).
  * [44] Shreejan Ram Shrestha, Nguyen Thi Kim Oanh, Quishi Xu, Maheswar Rupakheti, and Mark G Lawrence.  “Analysis of the vehicle fleet in the kathmandu valley for estimation of environment and climate co-benefits of technology intrusions” (2013).
  * [45] Guy Zisling, Lea Santos, and Yevgeny Bar Lev.  “How many particles make up a chaotic many-body quantum system?”.  SciPost Physics10 (2021).
  * [46] Christoph Berke, Evangelos Varvelis, Simon Trebst, Alexander Altland, and David P. DiVincenzo.  “Transmon platform for quantum computing challenged by chaotic fluctuations”.  Nature Communications13 (2022).
  * [47] Michael A Nielsen and Isaac L Chuang.  “Quantum computation and quantum information”.  Phys. Today 54, 60 (2001).
  * [48] Yilun Yang, Arthur Christianen, Mari Carmen Bañuls, Dominik S. Wild, and J. Ignacio Cirac.  “Phase-sensitive quantum measurement without controlled operations” (2023).  arXiv:2308.10796.
  * [49] S Lerma-Hernández, D Villaseñor, MA Bastarrachea-Magnani, EJ Torres-Herrera, Lea F Santos, and JG Hirsch.  “Dynamical signatures of quantum chaos and relaxation time scales in a spin-boson system”.  Physical Review E 100, 012218 (2019).
  * [50] Alaina M. Green, A. Elben, C. Huerta Alderete, Lata Kh Joshi, Nhung H. Nguyen, Torsten V. Zache, Yingyue Zhu, Bhuvanesh Sundar, and Norbert M. Linke.  “Experimental measurement of out-of-time-ordered correlators at finite temperature”.  Physical Review Letters128 (2022).