# Spectroscopic Confirmation of an Ultra-Massive Galaxy in a Protocluster at
z$\sim 4.9$∗

Stephanie M. Urbano Stawinski${}^{1,^{\dagger}}$ M. C. Cooper1 Ben Forrest2
Adam Muzzin3 Danilo Marchesini4 Gillian Wilson5 Percy Gomez6 Ian McConachie7
Z. Cemile Marsan3 Marianna Annunziatella8 Wenjun Chang7 1Department of Physics
& Astronomy, University of California, Irvine, 4129 Reines Hall, Irvine, CA
92697, USA  2Department of Physics and Astronomy, University of California
Davis, One Shields Avenue, Davis, CA 95616 USA  3Department of Physics and
Astronomy, York University, 4700, Keele Street, Toronto, ON MJ3 1P3, Canada
4Department of Physics and Astronomy, Tufts University, 574 Boston Avenue,
Medford, MA 02155, USA  5Department of Physics, University of California
Merced, 5200 North Lake Rd., Merced, CA 95343, USA  6W.M. Keck Observatory,
65-1120 Mamalahoa Hwy., Kamuela, HI 96743, USA  7Department of Physics &
Astronomy, University of California Riverside, 900 University Ave., Riverside,
CA 92521 USA  8Centro de Astrobiología (CSIC-INTA), Ctra de Torrejón a
Ajalvir, km 4, E-28850 Torrejón de Ardoz, Madrid, Spain

###### Abstract

We present spectroscopic confirmation of an ultra-massive galaxy (UMG) with
$\log(M_{\star}/M_{\odot})=10.98\pm 0.09$ at $z_{\mathrm{spec}}=4.8947$ in the
Extended Groth Strip (EGS), based on deep observations of Ly$\alpha$ emission
with Keck/DEIMOS. The ultra-massive galaxy (UMG-28740) is the most massive
member in one of the most significant overdensities in the EGS, with four
additional photometric members with $\log(M_{\star}/M_{\odot})>10.5$ within
$R_{\mathrm{proj}}\sim 1$ cMpc. Spectral energy distribution (SED) fitting
using a large suite of star formation histories and two sets of high-quality
photometry from ground- and space-based facilities consistently estimates the
mass of this object to be $\log(M_{\star}/M_{\odot})\sim 11$ with a small
standard deviation between measurements ($\sigma=0.09$). While the best-fit
SED models agree on stellar mass, we find discrepancies in the estimated star
formation rate for UMG-28740, resulting in either a star-forming or quiescent
system. JWST/NIRCam photometry of UMG-28740 strongly favors a quiescent
scenario, demonstrating the need for high-quality mid-IR observations.
Assuming the galaxy to be quiescent, UMG-28740 formed the bulk of its stars at
$z>10$ and is quenching at $z\sim 8$, resulting in a high star formation
efficiency at high redshift ($\epsilon\sim 0.2$ at $z\sim 5$ and
$\epsilon\gtrsim 1$ at $z\gtrsim 8$). As the most massive galaxy in its
protocluster environment, UMG-28740 is a unique example of the impossibly
early galaxy problem.

###### keywords:

galaxies:high-redshift – galaxies:evolution – galaxies:clusters

††thanks: †Corresponding author<EMAIL_ADDRESS>∗The data
presented herein were obtained at the W. M. Keck Observatory, which is
operated as a scientific partnership among the California Institute of
Technology, the University of California and the National Aeronautics and
Space Administration. The Observatory was made possible by the generous
financial support of the W. M. Keck Foundation.

##  1 Introduction

The discovery of extremely massive, passive galaxies at $z\sim 3-4$ has
challenged our understanding of galaxy formation and evolution at early times
(Glazebrook et al., 2017; Schreiber et al., 2018; Forrest et al., 2020a, b;
Valentino et al., 2020). These quiescent ultra-massive galaxies (UMGs,
$\hbox{${M}_{\star}$}>10^{11}~{}\hbox{${\rm M}_{\odot}$}$) have little to no
ongoing star formation, requiring a period of intense, bursty star formation
followed by rapid quenching within the first 2 Gyr of the Universe. Such
extreme systems are difficult to reproduce in galaxy formation models, with
current simulations struggling to form UMGs so early in cosmic time and and
often failing to quench these systems by $z\sim 3$ (Forrest et al., 2020a;
Valentino et al., 2023; Xie et al., 2024). Thus, although rare, confirmation
of UMGs provide powerful constraints on theoretical models of galaxy formation
and evolution (e.g. Baugh, 2006; Croton et al., 2006; Nelson et al., 2018).

Intriguingly, the modeled ages and star formation histories of UMGs at $z\sim
3$ suggest that some systems assembled the majority of their mass and
subsequently quenched as early as $z\sim 5$ (Forrest et al., 2020b; Schreiber
et al., 2018). Overall, there appears to be rapid evolution of the UMG
population at $4<z<6$, as some of these massive systems transition from highly
star-forming to passive. Imaging campaigns that originally uncovered UMGs at
$z\sim 3$ also corroborate these findings, yielding samples of photometric UMG
candidates at $z>4$ (e.g. Muzzin et al., 2013a, b; Weaver et al., 2023).
Meanwhile, measurements of the number density of UMGs find a significant
increase in their abundance with cosmic time from $z\sim 6$ to $z\sim 3$
(Stefanon et al., 2015; Davidzon et al., 2017; Marsan et al., 2022; Weibel et
al., 2024). Altogether, deep photometric surveys point to not only the
existence of these extreme systems at $z\sim 5$, but also that this epoch is
critical for understanding the rapid evolution of this extreme population.

While it is presumed that ultra-massive galaxies should exist at $z\sim 5$,
spectroscopic confirmation of candidate systems has been challenging. Studies
searching for UMGs at $z\sim 5$ have encountered difficulty in identifying and
detecting these systems via both ground- and space-based spectroscopy (Antwi-
Danso et al. 2023; Marsan et al. in prep). The few existing,
spectroscopically-confirmed UMGs at $z>4$ represent a biased population of
systems with exceptionally high star formation rates (${\rm SFR}\gtrsim
500~{}\hbox{${\rm M}_{\odot}$}~{}{\rm yr}^{-1}$) and with relatively poorly-
constrained stellar masses (e.g. Xiao et al., 2023, but see also Carnall et
al. 2023).

Herein, we present the discovery of a spectroscopically-confirmed, ultra-
massive ($\hbox{${M}_{\star}$}\sim 10^{11}~{}\hbox{${\rm M}_{\odot}$}$) galaxy
at $z=4.8947$ that resides within an extreme overdense region in the Extended
Groth Strip (EGS). In §2, we describe our target selection criterion and
spectroscopic observations. We present the resulting spectral fitting and
characterization of the protocluster environment in §3. We discuss the
competing star formation histories (SFH), compare our object to existing
simulations, and investigate the mass assembly history of this object in §4.
Finally, we summarize our results in §5. Unless otherwise specified, we adopt
a $\Lambda$CDM cosmology throughout with $\Omega_{\rm M}=0.3$,
$\Omega_{\Lambda}=0.70$, and $H_{0}=68$ km s-1 Mpc-1 (Planck Collaboration et
al., 2016).

##  2 Data

###  2.1 Target Selection and Photometric Measurements

We selected UMG candidates ($\hbox{${M}_{\star}$}\gtrsim 10^{11}~{}\hbox{${\rm
M}_{\odot}$}$) at $3.5<z_{\rm phot}<6.5$ in the Extended Groth Strip (EGS)
using the multi-wavelength photometric catalog from the Cosmic Assembly Near-
infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011;
Koekemoer et al. 2011). The CANDELS photometric catalog (Stefanon et al.,
2017), which spans $0.4-8\mu$m via a suite of ground- and space-based imaging,
incorporates photometric redshifts and stellar masses from seven separate
photo-$z$ codes and/or spectral energy distribution (SED) template sets. For
selection purposes, we used the median $z_{\mathrm{phot}}$ and ${M}_{\star}$
reported in the catalog. To increase the likelihood of detection of
Ly$\alpha$, we prioritized objects that had an apparent CFHT $r$-band AB
magnitude of $<28$ and star formation rate (SFR) of $\gtrsim 100~{}\hbox{${\rm
M}_{\odot}$}~{}{\rm yr}^{-1}$. To increase our sample size, we included UMG
candidates within the given redshift range that had lower SFRs and masses as
filler targets.

UMG-28740 ($\alpha=14^{\text{h}}19^{\text{m}}39.7^{\text{s}}$ and
$\delta=+52^{\circ}56^{\text{m}}56.5^{\text{s}}$), in particular, was selected
as a star-forming candidate, with an estimated SFR of $20-200~{}\hbox{${\rm
M}_{\odot}$}~{}{\rm yr}^{-1}$ and an apparent CFHT $r$-band magnitude of
$27.9^{+0.6}_{-0.4}$. Within the Stefanon et al. (2017) catalog, UMG-28740 has
a median photometric redshift of $z_{\mathrm{phot}}=4.495$ and a median
stellar mass of $\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})=11.0$. In
total, the CANDELS catalog reports photometric measurements of UMG-28740 at
$ugriz$ and $JHK_{s}$ bands from CFHT, F606W, F814W, F125W, F140W, and F160W
from HST, as well as all four Spitzer/IRAC channels.

In addition to the photometry from CANDELS, this galaxy was recently observed
with JWST/NIRCam (Rieke et al., 2003; Beichman et al., 2012) as part of the
Cosmic Evolution Early Release Survey (CEERS, PID 12117; Finkelstein et al.
2023; Bagley et al. 2023). Using JWST/NIRCam (pointing 2), UMG-28740 was
observed with the F115W, F150W, F200W, F277W, F356W, F410M, and F444W filters
at a $5\sigma$ depth ranging from AB magnitudes of $28.35-29.2$. While the
CEERS photometry was not incorporated in the initial target selection, it is
utilized in this work as a supplemental set of photometry in $\S\ref{sec:SED
Fit}$.

###  2.2 Spectroscopic Observations and Reduction

Using Keck/DEIMOS (Faber et al., 2003), we observed UMG-28740 in March and May
of 2023 on three different slitmasks. For all observations, we used the 600
lines mm-1 grating blazed at 7500 Å and tilted to a central wavelength of 7200
Å, with the GG455 order-blocking filter employed and no dithering between
exposures. With slit widths of $1^{\prime\prime}$, the spectral resolution
(FWHM) for the 600g grating on DEIMOS is $\sim 3.5$ Å (Weiner et al., 2006),
with a dispersion of $0.65$ Å per pixel. Calibrations for each slitmask
included three internal quartz lamp flat-field frames and an arc lamp spectrum
(using Kr, Ar, Ne, and Xe lamps). We were able to observe UMG-28740 with the
three different slitmasks for a combined exposure time of 6.9 hours. The
average seeing for these observations ranged from $\sim
0.8$-$1^{\prime\prime}$ with variable transparency.

Figure 1: (_top_) A $\sim 150$ Å-wide section of the co-added 2D spectrum for
UMG-28740 from our Keck/DEIMOS observations, highlighting the Lyman-$\alpha$
emission line. (_bottom_) The corresponding co-added 1D spectrum with the
resulting best-fit redshift shown in red.

The Keck/DEIMOS observations were reduced using the
spec2d111https://sites.uci.edu/spec2d/ DEEP2/DEEP3 DEIMOS data reduction
pipeline (Newman et al., 2013; Cooper et al., 2012), yielding sky-subtracted,
wavelength-calibrated 1D and 2D spectra for each target. We then combined the
three 1D spectra of UMG-28740 from the reduction of each slitmask, creating a
final co-added 1D spectrum as shown in Figure 1. For visualization of the
line, we also co-added the corresponding 2D spectra, as shown in the top panel
of Figure 1.

To measure a best-fit redshift, we utilized a custom template fitter (as used
in Urbano Stawinski et al. 2024) that incorporates both an emission-line
galaxy template as well as an asymmetric Gaussian profile to probe for a
single Ly$\alpha$ line. For UMG-28740, we found a single emission line nearby
the expected wavelength for Ly$\alpha$ given the photometric redshift. Fitting
this profile, we find the resulting heliocentric-corrected spectroscopic
redshift for UMG-28740 is $z_{\mathrm{spec}}=4.8947$ (best-fit model is shown
by the red line in Fig. 1).

While this redshift is determined by the presence of a single emission line,
there are multiple arguments consolidating this feature to be Ly$\alpha$.
First, the asymmetric profile of the emission line follows the characteristic
shape of Ly$\alpha$ emitters at $z>3$ (e.g. Verhamme et al., 2008; Martin et
al., 2015). Second, the measured redshift is consistent with two sets of
photometry from CANDELS and CEERS, including high-quality photometric
measurements from HST/WFC3 and JWST/NIRCam (see Fig. 2). Lastly, the location
of this galaxy coincides with a large overdense region at $z\sim 5$ (first
highlighted by Naidu et al. 2022 and later discussed in §3.2), which further
solidifies $z_{\mathrm{spec}}=4.8947$ as the correct redshift.

##  3 Analysis

###  3.1 SED Fitting

Figure 2: (_top_) From left to right: image stamps of UMG-28740 in the
JWST/NIRCam 1.50$\mu$m, 2.77$\mu$m, and 4.44$\mu$m filters, followed by an RGB
image (F150M in blue, F277M in green, and F444W in red). (_bottom_) Best-fit
SEDs (star-forming in blue and quenched in orange) from BAGPIPES forced at the
measured $z_{\mathrm{spec}}$. The two sets of observed photometry are shown as
dark blue (Stefanon et al., 2017) and orange (Bagley et al., 2023) markers.
The width of each SED indicates the $\pm 1\sigma$ confidence intervals
resulting from BAGPIPES. The open blue crosses show the predicted JWST/NIRCam
flux at F277W, F356W, F410M, and F444W for the blue star-forming SED. Here,
the predicted JWST/NIRCam flux for the star-forming scenario is under-
estimated at $2.77\mu{\rm m}$ and over-estimated at $3.56-4.44\mu{\rm m}$,
suggesting that UMG-28740 is likely quenched.

With our measured spectroscopic redshift ($z_{\mathrm{spec}}=4.8947$), we re-
fit the SED of UMG-28740. Here, we employ both the photometry from CANDELS
(Stefanon et al., 2017) and from CEERS (Bagley et al., 2023). From the CANDELS
photometric catalog, we exclude the HST/ACS photometry at $F606W$ and $F814W$
due to uncertainties related to sky subtraction, such that the SED fits are
made to the photometry at $ugrizJHK_{s}$ from CFHT, F125W, F140W, and F160W
from HST/WFC3-IR, and all four Spitzer/IRAC channels. Meanwhile, the CEERS
catalog includes photometric measurements at F606W and F814W from HST/ACS,
F125W, F140W, and F160W from HST/WFC3-IR, and finally F115W, F150W, F200W,
F277W, F356W, F410M, and F444W from JWST/NIRCam.

The majority of our SED analysis was completed using the SED-fitting code
Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation
or BAGPIPES222https://github.com/ACCarnall/bagpipes (Carnall et al., 2018).
BAGPIPES assumes Bruzal and Charlot stellar populations (Bruzual & Charlot,
2003), a Kroupa IMF (Kroupa & Boily, 2002), and a Calzetti dust law (Calzetti
et al., 2000). Adopting a double-power law star formation history, we measure
best-fit model SEDs for each set of photometry (from CANDELS and CEERS). The
SEDs are fit to 10 free parameters, allowing for variation in metallicity,
age, the slope of the Calzetti dust law, and the ionization parameter. The
full set of free parameters used and their allowed distributions are outlined
in Table 1.

Table 1: List of free parameters used in the SED fits with BAGPIPES Parameter | Unit | Range | Prior
---|---|---|---
Generala
Stellar Mass Formed | $\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$}$ | (1, 15) | Logarithmic
Peak Star Formation |  $\tau/$Gyr | (0.3, $t_{\mathrm{obs}}$) | Logarithmic
Age |  $t_{\mathrm{age}}/$Gyr | (0.1, 1) | Logarithmic
Metallicity | $Z/\hbox{${\rm Z}_{\odot}$}$ | (0.01, 2.5) | Logarithmic
$V$-band attenuation | $A_{V}$ | (0, 8) | Uniform
Deviation from | $\delta$ | (-0.3, 0.3) | Gaussianb
Calzetti slope
2175Å bump strength | $B$ | (0, 5) | Uniform
Ionization | $U$ | (-4, -2) | Logarithmic
Double-power lawc
double-power law | $\alpha$ | (0.01, 1000) | Logarithmic
falling slope
double-power law | $\beta$ | (0.01, 1000) | Logarithmic
rising slope
Burstd
Stellar Mass Formed | $\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$}$ | (0, 15) | Logarithmic
Age |  $t_{\mathrm{age}}/$Gyr | (0.1, 1) | Logarithmic
Metallicity | $Z/\hbox{${\rm Z}_{\odot}$}$ | (0.01, 2.5) | Logarithmic
AGNd
Continuum Flux at | $f_{5100}$ | (0, $10^{-19}$) | Uniform
5100Å | erg/s/cm2/Å
H$\alpha$ Flux | $f_{\mathrm{H\alpha}}$ | (0, $10^{-16.6}$) | Uniform
erg/s/cm2
H$\alpha$ Velocity | $\sigma_{\mathrm{H\alpha}}$/km/s | (1000, 5000) | Logarithmic
Dispersion
Power law slope: |  |  |
at $\lambda<5000$Å | $\alpha_{\lambda<5000}$ | (-2, 2) | Gaussiane
at $\lambda>5000$Å | $\alpha_{\lambda>5000}$ | (-2, 2) | Gaussianf
aParameters used in all SED fits.
bThe Gaussian distribution for the Calzetti slope has a mean of $\mu=0$
and a standard deviation of $\sigma=0.1$
cParameters specifically used in the double-power law fits.
dParameters used in added AGN and burst components.
eThe Gaussian distribution for the power law slope at at $\lambda<5000$Å
has a mean of $\mu=-1.5$ and a standard deviation of $\sigma=0.5$
fThe Gaussian distribution for the power law slope at at $\lambda>5000$Å
has a mean of $\mu=0.5$ and a standard deviation of $\sigma=0.5$

As shown in Figure 2, the resulting SED fits give a stellar mass of
$\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})=10.98^{+0.04}_{-0.03}$
using the CEERS photometry and $\log(\hbox{${M}_{\star}$}/\hbox{${\rm
M}_{\odot}$})=11.08^{+0.08}_{-0.1}$ using the CANDELS photometry. While the
stellar mass is tightly constrained between the two sets of photometry, we
find that the SFR is dependant on the photometric catalog utilized. The
CANDELS photometry yields a dusty, highly star-forming galaxy (${\rm
SFR}=285.0^{+125.5}_{-79.5}~{}\hbox{${\rm M}_{\odot}$}~{}{\rm yr}^{-1}$) and
the fit to the JWST/NIRCam photometry produces a very low star formation rate
(${\rm SFR}=9.5^{+4.8}_{-2.7}~{}\hbox{${\rm M}_{\odot}$}~{}{\rm yr}^{-1}$) and
low dust attenuation. The latter, nearly quenched SED is shown in orange in
Figure 2, and is consistent with both sets of photometry within $\pm 1\sigma$
error. On the other hand, the star-forming SED (shown in blue in Fig. 2) has
insufficient emission in the rest-frame optical ($\sim 4000$-$6000$Å), so as
to under-produce the flux measured in the JWST/NIRCam F277W band. In addition,
the star-forming SED overpredicts the observed flux in the JWST/NIRCam F356W,
F410M, and F444W bands due to contribution from strong emission lines such as
H$\alpha$ at $\sim 3.9\mu{\rm m}$. Based on the photometry alone, this favors
the quenched scenario for UMG-28740 (see §4.1 for a more detailed discussion
of these discrepancies).

To test the robustness of our measurements for mass and SFR, we ran BAGPIPES
again using a variety of different models; varying the included free
parameters, changing the assumed star formation history, restricting the
assumed metallicity to solar, and adding different components such as bursts
and an active galactic nucleus (AGN). In summary, we highlight the results
from five different models: a double-power law SFH, a double-power law SFH
including an AGN, a burst, or fixing the metallicity to solar metallicity, and
a delayed-$\tau$ model (see Table 2). Adding an AGN or burst component is of
interest to the evolution of massive galaxies in the early Universe, as the
spectra of a subset of UMGs have significant AGN contributions (Marsan et al.
2017; Forrest et al. 2020b, Forrest et al. in prep.) and current theoretical
models speculate that the rapid formation of such systems require extreme
bursts of star formation (Sun et al., 2023).

We also compared the results from BAGPIPES to fits using the C++-adapted SED
fitting code FAST++333https://github.com/cschreib/fastpp (Kriek et al., 2009).
With FAST++, we adopt the same parameters used in the literature to model UMGs
at $z\sim 3.5$ (Forrest et al., 2020a, b), assuming a Bruzal Charlot stellar
population (Bruzual & Charlot, 2003), the Calzetti dust law (Calzetti et al.,
2000), and a Chabrier IMF (Chabrier, 2003). The Chabrier IMF is related to the
Kroupa IMF used in BAGPIPES by a factor of $\sim 1$ (Wright et al., 2017),
allowing for direct comparison between results of the two SED codes. We allow
FAST++ to vary the extinction from $0<A_{V}<4$ with a $\Delta A_{V}=0.1$, the
age from $7<\log(\mathrm{Age})<$ age of the Universe at $z_{\mathrm{spec}}$
with $\Delta\log(\mathrm{Age})=0.05$, and a metallicity of $Z/\hbox{${\rm
Z}_{\odot}$}=0.02$ as in Forrest et al. (2020a, b). For the assumed SFH, we
adopt a double-exponential SFH as used in Forrest et al. (2020a, b).

Remarkably, regardless of the photometric dataset used, SFH assumed,
underlying SED code, or range of free parameters, we find all best-fit SEDs
result in a massive galaxy with a stellar mass of
$\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})\sim 11$. Across the
various fits, UMG-28740 has a median mass of
$\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})=10.98$ with a $1\sigma$
standard deviation of $0.09$ (see Table 2 for measurements from each SED fit).
Hereinafter, we conservatively use this stellar mass and error as our best
estimate for UMG-28740. This SED fitting analysis also highlights a
discrepancy between SFR and dust extinction – resulting in the two possible
scenarios, star-forming or quenched (see §4.1 for a more detailed discussion).

###  3.2 An Overdense Region around UMG 28740

UMG-28740 is embedded within a significant overdensity at $z\sim 5$, as traced
by the surface density of galaxies at $4.5<z_{\rm phot}<5.5$ in the CANDELS
photometric catalog. As shown in Figure 3, UMG-28740 resides near the center
of the overdense region and is the most massive member. However, we caution
that the observed overdensity is located near the edge of the CANDELS survey
area within the EGS, with UMG-28740 only $\sim 16^{\prime\prime}$ from the
edge of the HST/WFC3 imaging footprint, such that the full extent of the
structure may be underestimated. Naidu et al. (2022) first identified this
particular overdense region, finding that a $\sim 0.8~{}{\rm arcmin}^{2}$ area
centered on a lower-mass photometric neighbor of UMG-28740 contained $4\times$
the average surface density of sources at $4.5<z_{\rm phot}<5.5$. However,
this lower-mass neighbor ($\hbox{${M}_{\star}$}\sim 10^{9.6}~{}\hbox{${\rm
M}_{\odot}$}$, spectroscopically confirmed by Arrabal Haro et al. 2023b), is
located $\gtrsim 20^{\prime\prime}$ to the south ($\gtrsim 0.8~{}{\rm cMpc}$
at $z=5$) of UMG-28740 and outside of the apparent center of the overdense
region as traced by the CANDELS photometric catalog (see the highlighted red
cross in Fig. 3).

To investigate the extent of the overdense region around UMG-28740, we search
for other overdensities across the entire CANDELS footprint within the EGS
field by comparing the galaxy surface densities in redshift slices of $\Delta
z=1$ to the average surface density at each corresponding epoch. We use the
CANDELS catalog to identify sources at $2.5<z_{\mathrm{phot}}<6.5$, as this
catalog covers a larger footprint within the EGS (versus CEERS). We limit this
analysis to sources with $\log(\hbox{${M}_{\star}$}/\hbox{${\rm
M}_{\odot}$})>9$, so as to define a tracer population that is unbiased
(complete) over the full redshift range probed. Across the entire field, we
compute the galaxy surface density as a function of position on the sky within
an aperture of $45^{\prime\prime}\times 45^{\prime\prime}$ and a depth of
$\Delta z=1$. Stepping through redshift from $z=3$ to $z=6$ with a step size
of $\Delta z=0.1$, we identify regions that are overdense compared to the
average surface density at the corresponding redshift. We find that the region
hosting UMG-28740 is one of the most extreme overdense regions within the EGS
at $3<z<6$, with a surface density that is $\sim 10\times$ higher than the
average surface density of the field at $z=4.9$.

Figure 3: Overdensity hosting UMG-28740 (blue star) at $z\sim 5$. The
underlying 2D histogram and contours show the galaxy surface density for
objects in the CANDELS catalog at $4.5<z_{\mathrm{phot}}<5.5$ and
$\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})>9$. The beam size (with a
diameter of $4.5^{\prime\prime}$) used to convolve the density map is shown in
the upper right corner. Cyan circles denote photometric sources found in the
CANDELS (Stefanon et al., 2017) and CEERS (Bagley et al., 2023) catalogs,
where white circles are low-mass galaxies only found in the deeper CEERS
catalog. The size of the circle indicates the estimated stellar mass of the
object. Red crosses and the orange diamond are spectroscopically-confirmed
sources from Arrabal Haro et al. (2023a) and Urbano Stawinski et al. (2024),
respectively. The highlighted red cross indicates the location of the lower-
mass member from Naidu et al. (2022).

At $z\sim 5$, the effective radii of protoclusters are estimated to extend up
to $R_{e}=5-10$ cMpc (Chiang et al., 2013). UMG-28740 is positioned nearby the
apparent center of the overdense region, such that we can estimate
$R_{\mathrm{proj}}$ from UMG-28740 to trace the effective radius of the
protocluster. However, we caution this region is also near the edge of the EGS
field, limiting observed galaxies past $R_{\mathrm{proj}}\sim 0.6$ cMpc to the
north-west and potentially skewing the true center of the overdensity.
Nevertheless, within $R_{\mathrm{proj}}=3.44$ cMpc from UMG-28740, there are
26 photometric protocluster members at $4.5<z<5.5$ and with
$\hbox{${M}_{\star}$}>10^{9}~{}\hbox{${\rm M}_{\odot}$}$ in the CANDELS
photometric catalog along with 160 members identified in the CEERS catalog,
pushing down to lower stellar masses
($\hbox{${M}_{\star}$}>10^{8}~{}\hbox{${\rm M}_{\odot}$}$). Additionally, out
of all the photometric members there are $4$ massive neighbors with
$10.5<\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})<11$ within just
29.5′′, or $\sim 1.13$ cMpc, of UMG-28740. The closest of these massive
members is only separated from UMG-28740 by $0.29$ cMpc. This makes the
protocluster surrounding UMG-28740 the most compact overdensity of massive
$\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})>10.5$ galaxies within the
entire EGS at $z\sim 5$ (according to the CANDELS catalog), with only two
other compact regions each containing two similarly massive galaxies at a
similar spatial scale ($R_{\mathrm{proj}}<3$ cMpc). Moreover, there are no
other instances of $>3$ massive galaxies within $R_{\mathrm{proj}}<3$ cMpc at
this redshift in the CANDELS catalog.

Table 2: Stellar Mass and SFR from Alternative SED Fits SFH | Code | Mass | SFR
---|---|---|---
|  | $\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})$ |  ${\rm M}_{\odot}$/yr
CANDELS (Stefanon et al., 2017)
Double-power lawa | BAGPIPES | 11.08 ${}^{0.08}_{0.1}$ | 285.0 ${}^{125.5}_{79.5}$
Double-power law+AGN | BAGPIPES | 10.97 ${}^{0.04}_{0.04}$ | 8.6 ${}^{3.9}_{2.4}$
Double-power law+burst | BAGPIPES | 10.98 ${}^{0.08}_{0.09}$ | 419.6 ${}^{156.7}_{162.2}$
Double-power law at $Z_{\odot}$ | BAGPIPES | 11.18 ${}^{0.07}_{0.09}$ | 299.9 ${}^{99.5}_{60.1}$
Delayed-$\tau$ | BAGPIPES | 10.92 ${}^{0.11}_{0.13}$ | 312.1 ${}^{159.1}_{106.2}$
Double Exponential | FAST++ | 10.95${}^{0.00}_{0.11}$ | 24.5${}^{16.6}_{3.3}$
CEERS (Bagley et al., 2023)
Double-power lawa | BAGPIPES | 10.98 ${}^{0.04}_{0.03}$ | 9.5 ${}^{4.8}_{2.7}$
Double-power law+AGN | BAGPIPES | 10.96 ${}^{0.04}_{0.04}$ | 0.0 ${}^{0.0}_{0.0}$
Double-power law+burst | BAGPIPES | 11.17 ${}^{0.03}_{0.04}$ | 18.1 ${}^{6.2}_{4.7}$
Double-power law at $Z_{\odot}$ | BAGPIPES | 11.17 ${}^{0.03}_{0.03}$ | 21.0 ${}^{9.9}_{19.0}$
Delayed-$\tau$ | BAGPIPES | 11.05 ${}^{0.02}_{0.03}$ | 32.2 ${}^{10.8}_{5.7}$
aSED fit is shown in Fig. 2.

This protocluster also hosts 5 spectroscopically-confirmed galaxies from
Keck/DEIMOS (Urbano Stawinski et al., 2024) and JWST/NIRSpec (Arrabal Haro et
al., 2023a) spectroscopy within $R_{\mathrm{proj}}=3.44$ cMpc (see Fig. 3).
The measured mass of these spectroscopic members range from
$\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})=9.18-9.95$ and span 0.75
cMpc to 2.10 cMpc from UMG-28740. In summary, all members within
$R_{\mathrm{proj}}=4$ cMpc, both spectroscopically confirmed and all
photometric candidates, are estimated to be less massive than UMG-28740.

Figure 4: Star formation histories for the two competing scenarios, star-
forming (blue) and quenched (orange), based on the best-fit SEDs from BAGPIPES
and FAST++ (see Table 2). The solid lines, as well as the shaded region, show
the best-fit and $\pm 1\sigma$ double-power law SFH from BAGPIPES (shown in
Fig. 2). The formation time ($t_{\rm form}$) corresponds to the epoch at which
50% of the stellar mass is formed for the double-power law models. The higher-
quality JWST/NIRCam photometry favors the quenched scenario and suggests that
UMG-28740 formed much of its stellar mass as early as $z\gtrsim 10$.

##  4 Discussion

###  4.1 Investigation into Competing Star Formation Histories

As mentioned in §3.1, the two sets of photometry yield two different SFH
scenarios, one with a highly star-forming and dusty SED versus one with low
star formation and little dust extinction. The high-quality JWST/NIRCam
photometry favors a nearly quenched, SED. In contrast, the photometry from
CANDELS favors a star-forming SED, with the exception of a quenched solution
when assuming a double-exponential SFH with FAST++ (as described in §3.1). We
show a variety of these SFHs, highlighting the best fit SEDs from the final
BAGPIPES run, in Figure 4.

Looking at the SED fits alone, the $\pm 1\sigma$ spread of the orange,
quenched SED is consistent with both sets of photometry within the errors.
However, there is a clear offset between the JWST/NIRCam photometry and the
blue, star-forming SED, with the largest discrepancy at F277W and F410M. To
test the consistency between the NIRCam photometry and the star-forming SED,
we predict the flux that we would expect from a mock F277W, F356W, F410M, and
F444W observation of the star-forming SED fit, using the ‘predict’ function
within the BAGPIPES code (see the open blue crosses in Fig. 2). For F277W, the
expected flux would be elevated from the continuum due to the emission lines
nearby $2.77\mu{\rm m}$; however, we find this predicted flux is insufficient
to reproduce the observed F277W flux even within $\pm 1\sigma$. For F356W and
F410M, the predicted flux is elevated to a greater extent, due to
contributions from strong emission lines such as H$\alpha$ at $\sim 3.9\mu{\rm
m}$. However, this added emission over-predicts the flux measured by the
NIRCam photometry (also seen, to a lesser extent, in the F444W predicted
flux). Taken together, this suggests that the quenched scenario is more
consistent with the observed photometry from both catalogs, and thus more
likely to be the true SFH of UMG-28740.

One argument against a quenched SFH is the detection of Ly$\alpha$ emission.
However, this may be explained by either AGN contribution or remnant star
formation consistent with a post-starburst SED. Assuming no AGN contamination,
we again use the ‘predict’ function in BAGPIPES to estimate the equivalent
width (EW) of Ly$\alpha$ in the various SED fits. For the nearly quenched
scenario, the SFR from all models without an AGN component ranges from 9.5
$\hbox{${\rm M}_{\odot}$}~{}{\rm yr}^{-1}$ to 32.2 $\hbox{${\rm
M}_{\odot}$}~{}{\rm yr}^{-1}$. Measuring the EW of Ly$\alpha$ from the model
SEDs (using the $\pm 1\sigma$ spread in the quenched models) yields an
EW${}_{\mathrm{Ly\alpha}}=16.5-42.7$Å. We then estimate the EW of the observed
Ly$\alpha$ line from Keck/DEIMOS. We find the continuum level just redward of
the line profile and calculate the EW using a custom python code. Because the
continuum level is so low, we conservatively call this measurement an upper
limit on the true EW. Performing this analysis yields an observed
EW${}_{\mathrm{Ly\alpha}}>34.0$Å, within the range of the predicted values for
the post-starburst SED. Hence, we cannot rule out the nearly quenched scenario
by detection of Ly$\alpha$ emission alone.

We perform the same analysis on the SEDs with AGN contribution. With AGN
contribution we find the $\pm 1\sigma$ spread for the quenched SEDs yields
EW${}_{\mathrm{Ly\alpha}}=29.9-31.7$Å, which is below the estimated observed
Ly$\alpha$ EW. This would suggest that AGN contribution alone may not explain
the presence of Ly$\alpha$ emission seen in the Keck/DEIMOS spectroscopy, and
that likely there is some remnant star formation (even with the presence of
AGN contamination).

The best-fit SEDs for UMG-28740 yield two different SFH scenarios, a highly
star-forming galaxy or a nearly quiescent one. However, due to the flux
required by the JWST/NIRCam photometry in the mid-IR (F277W, F356W, F410M, and
F444W) and the consistency in the EW of the observed and predicted Ly$\alpha$
emission line, we argue that the quenched SED is likely the true SFH of this
object.

###  4.2 Comparison to Simulations

Figure 5: SFR versus stellar mass of UMG-28740 for the two scenarios, star-
forming (blue star) and quenched (orange star). The black circles show
galaxies from the Illustris TNG-300 simulation at $z=5$. The red dashed line
and white squares trace the median of the galaxies in the TNG sample, while
the red dotted lines show the $\pm 1\sigma$ spread. We also include other
spectroscopically-confirmed massive galaxies at this epoch from the literature
as purple circles (Xiao et al., 2023) and a peach cross (Carnall et al.,
2023).

The number density of UMGs at $z\sim 5$ is extremely low, especially for more
quiescent systems (Stefanon et al., 2015; Davidzon et al., 2017; Marsan et
al., 2022; Weaver et al., 2023). In addition, due to various observational
challenges associated with observing these systems, there are only a handful
of spectroscopically-confirmed UMGs at this epoch in the literature (e.g.
Antwi-Danso et al., 2023; Xiao et al., 2023; Carnall et al., 2023). To compare
the SFR of UMG-28740 to a large, unbiased sample of galaxies, we plot SFR
versus stellar mass for UMG-28740 with simulated galaxies. Many simulations
struggle to produce massive systems at high redshift, and most fail to quench
these objects at early epochs. However, the Illustris TNG-300 simulation
(Nelson et al., 2019) has been able to form massive systems and subsequently
suppress SFR to reproduce quiescent UMGs up to $z\sim 3.5$ (Forrest et al.,
2020a; Schreiber et al., 2018; Glazebrook et al., 2017). Hence, we compare
UMG-28740 to simulated galaxies in TNG-300 at $z=5$ (see Fig. 5).

While TNG-300 reproduces massive quiescent systems at $z\sim 3.5$, as well as
massive star-forming galaxies at earlier epochs (Forrest et al., 2020a),
TNG-300 fails to create massive and quiescent galaxies like UMG-28740
(assuming a quenched scenario) and one other spectroscopically-confirmed
quenched system in the literature (Carnall et al., 2023) at $z\sim 5$. The
measured SFR of UMG-28740 for the quenched scenario is $\sim 2$ dex below that
of any other simulated galaxy at $\log(\hbox{${M}_{\star}$}/\hbox{${\rm
M}_{\odot}$})\sim 11$. On the other hand, the TNG-300 simulation does
reproduce massive, star-forming galaxies at this epoch, including those from
Xiao et al. (2023) and the star-forming scenario of UMG-28740. The lack of
massive quenched objects is likely a result of one of two possibilities: the
TNG-300 simulation fails to catch the cessation of star formation at this
epoch or that the size of TNG-300 is too small ($V=302.6^{3}$ Mpc3) to capture
a system as rare as UMG-28740 and its protocluster environment.

Figure 6: Mass assembly history of UMG-28740 for the best-fit star-forming
models (blue lines) and quiescent models (orange lines). We include the
double-power law fits (solid line), the fixed $Z={\rm Z}_{\odot}$ double-power
law fits (dashed dotted lines), and the delayed-$\tau$ fits (dotted lines).
The $\pm 1\sigma$ uncertainties for the double-power law fits are shown as the
shaded blue and orange regions. The green tracks and squares show the assembly
histories for spectroscopically-confirmed UMGs at $z\sim 3$ from Forrest et
al. (2020b). Spectroscopically-confirmed massive galaxies at $z\sim 5$ from
Xiao et al. (2023) and (Carnall et al., 2023) are shown as purple circles and
a peach cross, respectively. Finally, we include a sample of $z>5$ DSFGs (grey
open hexagons) from the literature and a sample of photometric massive galaxy
candidates from Labbé et al. (2023) (cyan diamonds). The grey shaded regions
show the stellar mass prohibited by $\Lambda$CDM cosmology, following
$\hbox{${M}_{\star}$}^{\mathrm{max}}=\epsilon
f_{b}{M_{\mathrm{halo}}}^{\mathrm{max}}$ where $\epsilon$ is the star
formation efficiency and $f_{b}$ is the cosmic baryon fraction. We include
three different thresholds, by assuming the area of the EGS at a 20$\%$
efficiency and a 100$\%$ efficiency as well as the area of the full sky and
100$\%$ efficiency.

Recently, Xie et al. (2024) used the GAEA semi-analytic model (De Lucia et
al., 2024), applied to a larger box size than that of TNG-300 ($V=685^{3}$
Mpc3), to study the formation of massive quiescent galaxies at early cosmic
times. In summary, Xie et al. (2024) find a larger relative number of quenched
massive galaxies at $z\sim 4.5$, with a quenched fraction for systems with
$\hbox{${M}_{\star}$}=(0.8-1.2)\times 10^{11}~{}\hbox{${\rm M}_{\odot}$}$ of
$f_{q}\sim 0.06$. Xie et al. (2024) suggest that the primary quenching
mechanism for massive galaxies at high redshift is accretion disk feedback via
quasar winds. Yet even within the GAEA simulation UMG-28740 is exceptionally
rare; galaxies as quiescent as the best-fit quenched model, with a specific
SFR of $\sim 10^{-10}~{}{\rm yr}^{-1}$, have a number density of only $\sim
2\times 10^{-7}$ per ${\rm cMpc}^{3}$ at $z\sim 5$ in GAEA (De Lucia et al.,
2024).

###  4.3 Mass Assembly History of UMG-28740

To investigate the potential evolution of UMG-28740, in Figure 6 we plot the
mass assembly history corresponding to the best-fit star-forming and quiescent
models from BAGPIPES (see Fig. 4). We compare the mass growth of UMG-28740 to
that of spectroscopically-confirmed UMGs at $z\sim 3.5$ from the Massive
Ancient Galaxies At $z>3$ NEar-infrared (MAGAZ3NE) survey (Forrest et al.,
2020a, b, 2022), a Keck/MOSFIRE effort that has successfully spectroscopically
confirmed a large sample of star-forming and quenched UMGs at $z\sim 3.5$. We
also include a sample of massive dusty star-forming galaxies (DSFGs) at $z>5$
from the literature (Riechers et al., 2020; Williams et al., 2019; Jin et al.,
2019; Zavala et al., 2018; Pavesi et al., 2018; Riechers et al., 2017;
Strandet et al., 2017; Ma et al., 2015; Cooray et al., 2014; Riechers et al.,
2013; Walter et al., 2012), including the dusty star-forming UMGs
spectroscopically confirmed by Xiao et al. (2023). Finally, we also include
the massive, spectroscopically-confirmed quiescent galaxy from Carnall et al.
(2023) and the photometric massive galaxy candidates from Labbé et al. (2023).
Below we discuss the possible connections between UMG-28740 and these various
populations.

For the mass assembly history predicted by the best-fit star-forming models,
UMG-28740 is consistent with being the progenitor of the existing UMG
population at $z\sim 3.5$, such that the ongoing star formation in the system
will continue to grow its stellar mass at $3\lesssim z\lesssim 5$, likely
becoming similar to one of the more massive UMGs detected by the MAGAZ3NE
survey. Similarly, the DSFG population at $z\sim 5$, including the galaxies
identified by Xiao et al. (2023), are also likely to evolve into members of
the existing UMG population at $z\sim 3.5$ (Forrest et al., 2020b). Looking
towards higher $z$ for the best-fit star-forming models, the assembly history
of UMG-28740 is consistent with the massive galaxy candidates from Labbé et
al. (2023) as potential progenitors, with UMG-28740 forming much of its
stellar mass at $z\lesssim 8$.

On the other hand, assuming the more-likely quenched SFH, UMG-28740 will
evolve into a lower-mass UMG at $z\sim 3.5$ ($\hbox{${M}_{\star}$}\sim
10^{11}~{}\hbox{${\rm M}_{\odot}$}$). Interestingly, however, in this scenario
UMG-28740 would not be included in the sample of known UMGs at $z\sim 3.5$.
Due to the difficulty in detecting absorption lines in faint systems from the
ground, MAGAZ3NE observations were limited to an AB magnitude of $K_{s}<21.7$.
In contrast, if we passively evolve the quenched model SED for UMG-28740 at
$z=4.8947$ to $z=3.5$, we find that UMG-28740 would have an AB magnitude of
$K_{s}=22.7$, fainter than the selection limit for the MAGAZ3NE survey and
indicating that existing samples of UMGs at intermediate redshift may be
biased against the oldest systems (see Forrest et al. 2020b). Instead,
UMG-28740 would more likely resemble much fainter UMGs at $z\sim 3$, such as
those that have been targeted with JWST/NIRSpec (Nanayakkara et al., 2024;
Glazebrook et al., 2024).

Towards high-$z$, the evolutionary history of UMG-28740 drastically changes
assuming a quenched model. At $z\sim 5$, UMG-28740 is already an extreme
system, more massive than the only other spectroscopically-confirmed quiescent
system at this epoch (Carnall et al., 2023). Moreover, in contrast to the
massive quiescent galaxy from Carnall et al. (2023), which formed its mass at
$z\sim 7$, we find UMG-28740 to already be quenching by $z\sim 8$, with the
bulk of its stellar mass formed at much earlier cosmic time ($z>10$). Based on
the quenched SED fits, the mass assembly of UMG-28740 is inconsistent with the
DSFG population at $z>5$ and the massive galaxy candidates from Labbé et al.
(2023), instead forming a more extreme amount of mass at very high redshift.

The formation of such massive galaxies at high redshift is a highly active
area of study thanks to a wealth of new observations from JWST. For example,
early JWST imaging campaigns discovered a surplus of massive galaxy candidates
at $z=8-14$ that stand in tension with predictions from $\Lambda$CDM (Casey et
al., 2023; Labbé et al., 2023; Boylan-Kolchin, 2023). These results have
sparked a myriad of questions and scrutiny in our current understanding of the
physical processes driving galaxy evolution, including those that underpin SED
fitting codes. The extreme build-up of mass required to assemble massive
galaxies suggests intense bursts of star formation that convert $>20\%$ of the
available baryons into stars (Steinhardt et al., 2016; Boylan-Kolchin, 2023;
Sun et al., 2023; Casey et al., 2023; Xiao et al., 2023; Antwi-Danso et al.,
2023).

To investigate the star formation efficiency of UMG-28740, we plot the
limiting mass as imposed by a $\Lambda$CDM cosmology (shown as the grey shaded
regions in Fig. 6). To do this, we follow the methodology from Boylan-Kolchin
(2023), first calculating the number density of halos of a given mass as a
function of redshift using hmf444https://hmf.readthedocs.io, a python package
that calculates the halo mass function at different redshifts (Murray, 2014).
We then find the maximum halo mass, given a survey volume, for which we would
observe $1$ such system. To convert to stellar mass, we follow
$\hbox{${M}_{\star}$}^{\mathrm{max}}=\epsilon
f_{b}{M_{\mathrm{halo}}}^{\mathrm{max}}$, where $f_{b}$ is the cosmic baryon
fraction and $\epsilon$ is the efficiency of converting baryons into stars. We
perform this calculation for the area of the CANDELS footprint within the EGS
(205 arcmin2, Stefanon et al. 2017) using a line-of-sight depth of $\Delta
z=2$ and an efficiency of 20$\%$ ($\epsilon=0.2$) and 100$\%$
($\epsilon=1.0$). We also calculate the maximum stellar mass for observations
across the full sky with an efficiency of 100$\%$, which is the upper limit
for what is allowed in a $\Lambda$CDM cosmology.

At $z\sim 5$, UMG-28740 has a stellar mass consistent with a relatively high
star formation efficiency ($\sim 20\%$), comparable to that found in larger
samples of massive galaxies at high redshift (e.g. Chworowsky et al., 2023).
In contrast, recent analyses of massive star-forming galaxies at $z>5$ find
even more extreme star formation efficiencies ($\epsilon\sim 0.5$, Xiao et
al., 2023; Casey et al., 2023). If UMG-28740 is quiescent at $z\sim 5$, then
its evolution is also consistent with an elevated formation efficiency at
higher $z$. As shown in Figure 6, at $z\gtrsim 8$ where the quenched model for
UMG-28740 is still star forming, the implied formation efficiency is on the
order of (and potentially in excess of) unity, consistent with the
expectations from models of feedback-free star formation (Dekel et al., 2023).
Moreover, for the quenched model, UMG-28740 reaches significant tension with
$\Lambda$CDM at even earlier times, exceeding $\epsilon=1$ across the area of
the full sky at $z\sim 15$. This suggests UMG-28740 formed its stars at nearly
impossible rates as currently predicted by theoretical models.

###  4.4 RUBIES-EGS-QG-1

While we finalized this manuscript, de Graaff et al. (2024) presented a
JWST/NIRSpec PRISM and G395M spectrum of UMG-28740 (identified as RUBIES-EGS-
QG-1 in de Graaff et al. 2024). The spectroscopic observations were collected
in March 2024, with coverage of the Balmer break and the H$\alpha$ emission
line. This effort was part of the RUBIES survey (GO-4233; PIs A. de Graaff and
G. Brammer), which uses JWST/NIRSpec to observe galaxies in the EGS and UDS
fields. In this section, we discuss the results of this independent work and
compare those to our work.

The physical properties of RUBIES-EGS-QG-1 were measured from SED fitting of
the JWST/NIRSpec PRISM spectrum, as well as photometry from HST/WFC3 (F125W,
F140W, and F160W), HST/ACS (F814W and F606W), and JWST/NIRCam (F115W, F150W,
F200W, F277W, F356W, F444W, and F410M). The spectrum and photometry were fit
simultaneously via Prospector (Leja et al., 2017), assuming a Chabrier IMF
(Chabrier, 2003), stellar population models from Flexible Stellar Population
Synthesis (FSPS, Conroy & Gunn, 2010; Conroy et al., 2009), the MILES spectral
library (Sánchez-Blázquez et al., 2006), and MESA Isochrones and stellar
tracks (MIST, Choi et al., 2016; Dotter, 2016). With Prospector, de Graaff et
al. (2024) explored four different star formation history models – a
delayed-$\tau$, rising, and fixed $Z=Z_{\odot}$ model as well as a SFH
comprised of 14 bins of constant star formation. The resulting best-fit
stellar mass is in agreement with our results, yielding
$\hbox{${M}_{\star}$}=8.7^{+0.5}_{-0.4}\times 10^{10}~{}\hbox{${\rm
M}_{\odot}$}$. With inclusion of mid-IR spectroscopy from JWST/NIRCam, de
Graaff et al. (2024) pin down the SFR in the last 100 Myr as ${\rm
SFR}_{100}=1.9^{+1.5}_{-1.0}~{}\hbox{${\rm M}_{\odot}$}~{}{\rm yr}^{-1}$. This
independent work uses a separate SED modeling code and includes a high-quality
spectrum that probes key spectral features, yet the final results are in
excellent agreement with our quenched model.

The complementary analysis from de Graaff et al. (2024), taken together with
the conclusions of our SED fitting, highlights the importance of high-quality
mid-IR photometry – currently only accessible via JWST/NIRCam – in
understanding systems like UMG-28740. Only with such data can we accurately
measure key physical properties of $z\gtrsim 4$ galaxies. This example
(UMG-28740) demonstrates the importance of JWST as we continue to search and
study the rarest and most massive systems in the Universe.

##  5 Conclusions

In this work, we present spectroscopic confirmation of an ultra-massive galaxy
(UMG, ID 28740) with $\log(\hbox{${M}_{\star}$}/\hbox{${\rm
M}_{\odot}$})=10.98\pm 0.09$ at $z_{\mathrm{spec}}=4.8947$ based on deep
spectroscopy with Keck/DEIMOS. UMG-28740 is one of the most massive
spectroscopically-confirmed galaxies at $z\sim 5$ and lies within one of the
most extreme overdensities in the EGS. Here, we summarize the main results of
our analysis:

  * •

Using the spectroscopic redshift, we run a large suite of SED modeling
(assuming various SFHs and using two SED fitting codes, BAGPIPES and FAST++)
on two sets of photometric data from CANDELS (Stefanon et al., 2017) and CEERS
(Bagley et al., 2023). Regardless of the SED code, set of free parameters, or
SFH assumed, the stellar mass of UMG-28740 is tightly constrained, with a
median mass of $\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})=10.98$ and
a standard deviation of $\pm 0.09$ between all runs.

  * •

We find a discrepancy between the measured SFR and photometric dataset
assumed, with multiple SED fits to the CANDELS photometry yielding a highly
star-forming object. However, when fitting to the JWST/NIRCam photometry from
CEERS, UMG-28740 is nearly quiescent, with SFR $<32.2~{}\hbox{${\rm
M}_{\odot}$}~{}{\rm yr}^{-1}$. Futhermore, while the quenched SED is
consistent with both sets of photometry, the predicted flux from the best-fit
star-forming SED is under-predicts the measured flux at F277W and over-
predicts the measured flux at F356W, F410M, and F444W; suggesting a quenched
scenario is more likely. Finally, we show the presence of Ly$\alpha$ does not
rule out a post-starburst, nearly quenched scenario.

  * •

We find UMG-28740 to be the most massive member of a protocluster environment
at $z\sim 5$, with the surface density of galaxies in the region around
UMG-28740 roughly $10\times$ higher than the average coeval field surface
density. Within $R_{\mathrm{proj}}=3.44$ cMpc of UMG-28740, we find 26
photometric members with $\hbox{${M}_{\star}$}>10^{9}~{}\hbox{${\rm
M}_{\odot}$}$ and 160 photometric members with
$\hbox{${M}_{\star}$}>10^{8}~{}\hbox{${\rm M}_{\odot}$}$ at $4.5<z<5.5$, using
the CANDELS and CEERS photometric catalogs, respectively. Finally, we find
four additional massive galaxies
($\hbox{${M}_{\star}$}>10^{10.5}~{}\hbox{${\rm M}_{\odot}$}$) within
$0.29-1.13$ cMpc of UMG-28740 in the CANDELS photometric catalog.

  * •

Comparing the SFR of UMG-28740 to the Illustris TNG-300 simulation box at
$z=5$ shows that, if quenched, UMG-28740 is extremely rare – with a SFR $\sim
2$ dex lower than any simulated galaxy at
$\log(\hbox{${M}_{\star}$}/\hbox{${\rm M}_{\odot}$})\gtrsim 11$.

  * •

We find UMG-28740 has a high star formation efficiency ($\sim 20\%$) at $z\sim
5$. Assuming a quenched model, this efficiency exceeds $\epsilon=1$ for the
area of the full sky, suggesting UMG-28740 would need to form its stars at
faster rates than predicted by theoretical models.

## Acknowledgements

Based on observations with the NASA/ESA/CSA James Webb Space Telescope
obtained at the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Incorporated, under
NASA contract NAS5-03127. The authors wish to recognize and acknowledge the
very significant cultural role and reverence that the summit of Maunakea has
always had within the indigenous Hawaiian community. We are most fortunate to
have the opportunity to conduct observations from this mountain. This work is
based in part on observations made with the NASA/ESA Hubble Space Telescope
and with the NASA/ESA/CSA James Webb Space Telescope from the Space Telescope
Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Inc., under NASA contracts NAS 5–26555 and NAS 5-03127.
SMUS and MCC acknowledge support from the National Science Foundation through
grant AST-1815475. BF acknowledges support from JWST-GO-02913.001-A. GW
gratefully acknowledges support from the National Science Foundation through
grant AST-2205189 and from HST program number GO-16300. This research was
supported by the International Space Science Institute (ISSI) in Bern, through
ISSI International Team project “Understanding the evolution and transitioning
of distant proto-clusters into clusters". We also gratefully acknowledge the
Lorentz Center in Leiden (NL) for facilitating discussions related to this
project. MCC thanks Michael Boylan-Kolchin, Kristin Turney, and Lucas Cooper
for helpful discussions related to this work. SMUS thanks Adam Carnall for
additional support with BAGPIPES functionality.

## References

  * Antwi-Danso et al. (2023) Antwi-Danso J., et al., 2023, arXiv e-prints, p. arXiv:2307.09590
  * Arrabal Haro et al. (2023a) Arrabal Haro P., et al., 2023a, Nature, 622, 707
  * Arrabal Haro et al. (2023b) Arrabal Haro P., et al., 2023b, ApJ, 951, L22
  * Bagley et al. (2023) Bagley M. B., et al., 2023, ApJ, 946, L12
  * Baugh (2006) Baugh C. M., 2006, Reports on Progress in Physics, 69, 3101
  * Beichman et al. (2012) Beichman C. A., Rieke M., Eisenstein D., Greene T. P., Krist J., McCarthy D., Meyer M., Stansberry J., 2012, in Clampin M. C., Fazio G. G., MacEwen H. A., Oschmann Jacobus M. J., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave. p. 84422N, doi:10.1117/12.925447
  * Boylan-Kolchin (2023) Boylan-Kolchin M., 2023, Nature Astronomy, 7, 731
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, ApJ, 533, 682
  * Carnall et al. (2018) Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379
  * Carnall et al. (2023) Carnall A. C., et al., 2023, Nature, 619, 716
  * Casey et al. (2023) Casey C. M., et al., 2023, arXiv e-prints, p. arXiv:2308.10932
  * Chabrier (2003) Chabrier G., 2003, Publications of the Astronomical Society of the Pacific, 115, 763
  * Chiang et al. (2013) Chiang Y.-K., Overzier R., Gebhardt K., 2013, ApJ, 779, 127
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, The Astrophysical Journal, 823, 102
  * Chworowsky et al. (2023) Chworowsky K., et al., 2023, arXiv e-prints, p. arXiv:2311.14804
  * Conroy & Gunn (2010) Conroy C., Gunn J. E., 2010, ApJ, 712, 833
  * Conroy et al. (2009) Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486
  * Cooper et al. (2012) Cooper M. C., Newman J. A., Davis M., Finkbeiner D. P., Gerke B. F., 2012, spec2d: DEEP2 DEIMOS Spectral Pipeline, Astrophysics Source Code Library, record ascl:1203.003 (ascl:1203.003)
  * Cooray et al. (2014) Cooray A., et al., 2014, ApJ, 790, 40
  * Croton et al. (2006) Croton D. J., et al., 2006, MNRAS, 365, 11
  * Davidzon et al. (2017) Davidzon I., et al., 2017, A&A, 605, A70
  * De Lucia et al. (2024) De Lucia G., Fontanot F., Xie L., Hirschmann M., 2024, arXiv e-prints, p. arXiv:2401.06211
  * Dekel et al. (2023) Dekel A., Sarkar K. C., Birnboim Y., Mandelker N., Li Z., 2023, MNRAS, 523, 3201
  * Dotter (2016) Dotter A., 2016, The Astrophysical Journal Supplement Series, 222, 8
  * Faber et al. (2003) Faber S. M., et al., 2003, in Iye M., Moorwood A. F. M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes. pp 1657–1669, doi:10.1117/12.460346
  * Finkelstein et al. (2023) Finkelstein S. L., et al., 2023, ApJ, 946, L13
  * Forrest et al. (2020a) Forrest B., et al., 2020a, ApJ, 890, L1
  * Forrest et al. (2020b) Forrest B., et al., 2020b, ApJ, 903, 47
  * Forrest et al. (2022) Forrest B., et al., 2022, ApJ, 938, 109
  * Glazebrook et al. (2017) Glazebrook K., et al., 2017, Nature, 544, 71
  * Glazebrook et al. (2024) Glazebrook K., et al., 2024, Nature, 628, 277
  * Grogin et al. (2011) Grogin N. A., et al., 2011, ApJS, 197, 35
  * Jin et al. (2019) Jin S., et al., 2019, ApJ, 887, 144
  * Koekemoer et al. (2011) Koekemoer A. M., et al., 2011, ApJS, 197, 36
  * Kriek et al. (2009) Kriek M., van Dokkum P. G., Labbé I., Franx M., Illingworth G. D., Marchesini D., Quadri R. F., 2009, ApJ, 700, 221
  * Kroupa & Boily (2002) Kroupa P., Boily C. M., 2002, MNRAS, 336, 1188
  * Labbé et al. (2023) Labbé I., et al., 2023, Nature, 616, 266
  * Leja et al. (2017) Leja J., Johnson B. D., Conroy C., van Dokkum P. G., Byler N., 2017, ApJ, 837, 170
  * Ma et al. (2015) Ma J., et al., 2015, ApJ, 812, 88
  * Marsan et al. (2017) Marsan Z. C., Marchesini D., Brammer G. B., Geier S., Kado-Fong E., Labbé I., Muzzin A., Stefanon M., 2017, The Astrophysical Journal, 842, 21
  * Marsan et al. (2022) Marsan Z. C., et al., 2022, The Astrophysical Journal, 924, 25
  * Martin et al. (2015) Martin C. L., Dijkstra M., Henry A., Soto K. T., Danforth C. W., Wong J., 2015, ApJ, 803, 6
  * Murray (2014) Murray S., 2014, HMF: Halo Mass Function calculator, Astrophysics Source Code Library, record ascl:1412.006
  * Muzzin et al. (2013a) Muzzin A., et al., 2013a, ApJS, 206, 8
  * Muzzin et al. (2013b) Muzzin A., et al., 2013b, ApJ, 777, 18
  * Naidu et al. (2022) Naidu R. P., et al., 2022, arXiv e-prints, p. arXiv:2208.02794
  * Nanayakkara et al. (2024) Nanayakkara T., et al., 2024, Scientific Reports, 14, 3724
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 475, 624
  * Nelson et al. (2019) Nelson D., et al., 2019, Computational Astrophysics and Cosmology, 6, 2
  * Newman et al. (2013) Newman J. A., et al., 2013, ApJS, 208, 5
  * Pavesi et al. (2018) Pavesi R., et al., 2018, ApJ, 861, 43
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Riechers et al. (2013) Riechers D. A., et al., 2013, Nature, 496, 329
  * Riechers et al. (2017) Riechers D. A., et al., 2017, ApJ, 850, 1
  * Riechers et al. (2020) Riechers D. A., et al., 2020, ApJ, 895, 81
  * Rieke et al. (2003) Rieke M. J., et al., 2003, in Mather J. C., ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4850, IR Space Telescopes and Instruments. pp 478–485, doi:10.1117/12.489103
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez P., et al., 2006, MNRAS, 371, 703
  * Schreiber et al. (2018) Schreiber C., et al., 2018, A&A, 618, A85
  * Stefanon et al. (2015) Stefanon M., et al., 2015, ApJ, 803, 11
  * Stefanon et al. (2017) Stefanon M., et al., 2017, ApJS, 229, 32
  * Steinhardt et al. (2016) Steinhardt C. L., Capak P., Masters D., Speagle J. S., 2016, ApJ, 824, 21
  * Strandet et al. (2017) Strandet M. L., et al., 2017, ApJ, 842, L15
  * Sun et al. (2023) Sun G., Faucher-Giguère C.-A., Hayward C. C., Shen X., Wetzel A., Cochrane R. K., 2023, The Astrophysical Journal Letters, 955, L35
  * Urbano Stawinski et al. (2024) Urbano Stawinski S. M., et al., 2024, MNRAS, 528, 5624
  * Valentino et al. (2020) Valentino F., et al., 2020, ApJ, 889, 93
  * Valentino et al. (2023) Valentino F., et al., 2023, ApJ, 947, 20
  * Verhamme et al. (2008) Verhamme A., Schaerer D., Atek H., Tapken C., 2008, A&A, 491, 89
  * Walter et al. (2012) Walter F., et al., 2012, Nature, 486, 233
  * Weaver et al. (2023) Weaver J. R., et al., 2023, A&A, 677, A184
  * Weibel et al. (2024) Weibel A., et al., 2024, arXiv e-prints, p. arXiv:2403.08872
  * Weiner et al. (2006) Weiner B. J., et al., 2006, ApJ, 653, 1027
  * Williams et al. (2019) Williams C. C., et al., 2019, ApJ, 884, 154
  * Wright et al. (2017) Wright A. H., et al., 2017, Monthly Notices of the Royal Astronomical Society, 470, 283
  * Xiao et al. (2023) Xiao M., et al., 2023, arXiv e-prints, p. arXiv:2309.02492
  * Xie et al. (2024) Xie L., et al., 2024, arXiv e-prints, p. arXiv:2402.01314
  * Zavala et al. (2018) Zavala J. A., et al., 2018, Nature Astronomy, 2, 56
  * de Graaff et al. (2024) de Graaff A., et al., 2024, arXiv e-prints, p. arXiv:2404.05683


# Benchmarking Quantum Annealers with Near-Optimal Minor-Embedded Instances

Valentin Gilbert  Université Paris-Saclay
CEA, List, F-91120
Palaiseau, France
<EMAIL_ADDRESS>Julien Rodriguez  Université de Montpellier
LIRMM, CNRS
Montpellier, France
<EMAIL_ADDRESS>Stéphane Louise  Université Paris-Saclay
CEA, List, F-91120
Palaiseau, France
<EMAIL_ADDRESS>
###### Abstract

Benchmarking Quantum Process Units (QPU) at an application level usually
requires considering the whole programming stack of the quantum computer. One
critical task is the minor-embedding (resp. transpilation) step, which
involves space-time overheads for annealing-based (resp. gate-based) quantum
computers. This paper establishes a new protocol to generate graph instances
with their associated near-optimal minor-embedding mappings to D-Wave Quantum
Annealers (QA). This set of favorable mappings is used to generate a wide
diversity of optimization problem instances. We use this method to benchmark
QA on large instances of unconstrained and constrained optimization problems
and compare the performance of the QPU with efficient classical solvers. The
benchmark aims to evaluate and quantify the key characteristics of instances
that could benefit from the use of a quantum computer. In this context,
existing QA seem best suited for unconstrained problems on instances with
densities less than $10\%$. For constrained problems, the penalty terms used
to encode the hard constraints restrict the performance of QA and suggest that
these QPU will be less efficient on these problems of comparable size.

###### Index Terms:

Benchmark, Quantum Annealer, Optimization, QUBO, Ising model, maxcut, maximum
independent set

##  I Introduction

The past five years have seen the rise of quantum computers with successfully
implemented quantum supremacy experiments. These experiments were demonstrated
with superconducting qubits [1] and with photonic qubits [2, 3]. More
recently, D-Wave claimed that their QA reached quantum supremacy [4]. These
demonstrations constitute strong experimental steps towards achieving a useful
quantum advantage but are insufficient to demonstrate the utility of quantum
computers to solve real-world and useful problems. Combinatorial optimization
problems constitute a set of relevant problems to evaluate the potential
quantum advantage that could be brought by Noisy Intermediate Scale Quantum
(NISQ) process units. The Q-Score is an example of application-oriented
protocol designed to benchmark QPU using random max-cut instances based on
Erdős-Rényi graphs[5]. Several frameworks have been recently designed to
benchmark annealing-based and gate-based models in a holistic way [6]. One
example is the QUARK framework, which explores real industrial use cases [7].
Various national and international initiatives are developing application-
oriented frameworks, including the QED Consortium [8] and the BACQ project[9].
These frameworks are completed with attempts to quantify the coverage of
instance sets to gauge the fairness of a benchmark experiment [10, 11]. Recent
studies gather specific instance sets designed to benchmark quantum computers,
such as HamLib [12], which proposes various encodings of Hamiltonians for
optimization problems and condensed matter models. The MQTBench library [13]
establishes a database of quantum circuits that can be used for gate-based
QPU. One can also rely on previously developed reference libraries of problems
used for classical computers, such as the MQLib [14] and the QPLib [15]. The
interested reader can refer to a recent survey done on this subject [16].

The methods described above are based on the same execution scheme: an input
instance is defined and a preprocessing routine is used to map the problem on
the QPU. This task is called the compilation or transpilation for gate-based
quantum computers. It converts the original gate set to a physically
realizable gate set and adds SWAP gates to compensate for the limited
interconnects of the quantum chip. For QA such as D-Wave Systems [17], the
preprocessing step consists of finding a transformation that maps each
problem’s variable to a set of physical variables that can be straightly
encoded to the physical qubits of the quantum annealer. This task is called
minor-embedding [18]. Transpilation as well as minor-embedding tasks are hard
for classical computers and directly impact the quality of the solutions
provided by the QPU.

A second approach uses sets of crafted instances to benchmark the QPU [16].
These instances are designed to reduce the impact of the preprocessing step
and constitute friendly use cases for the QPU but remain hard for classical
optimizers. Table I shows recent experiments led on D-Wave QA using this type
of instance. Each experiment shows that QA can excel at solving specific sets
of instances compared to classical solvers. These instances have two common
features: they have many variables and are almost always subgraphs of the
qubit layout, which avoids using extra qubits to embed the problem on the
quantum chip. Crafted instances have been used to compare the performance of
IBM gate-based QPU to D-Wave QA, leading to the superiority of QA for solving
optimization problems [19].

TABLE I: Related work using crafted instance sets to demonstrate the performance of D-Wave quantum annealers. Source | Instance type | Instance topology | #Variables | #Qubits | Tested QPU | Minor-Embedding | Classical solvers
---|---|---|---|---|---|---|---
[20] (2021) | BFM∗, FBFM∗, CBFM∗ | Chimera graph | 2032 | 2032 | 2000Q | QPU chip sub-graph | Global & local search
[21] (2022) | CBFM∗-P | Pegasus graph | 5387 | 5387 | Advantage_System4.1 | QPU chip sub-graph | Global & local search
[22] (2023) | Unweighted Max-cut | 3-regular graph | 4-320 | 4-320 | Advantage_System4.1 | QPU chip sub-graph | Exact solver
[4] (2024) |  2D & 3D Ising spin glass |  square, cubic, diamond, biclique | 16-567 | 16-576 |  Advantage_System4.1 and Advantage2 |  QPU chip sub-graph 2 qubits/var (biclique) | Exact solver
∗BFM, FBFM and CBFM respectively relates to Biased Ferromagnet, Frustrated
Biased Ferromagnet and Corrupted Biased Ferromagnet Ising models.

This paper provides a new method to create large sets of crafted instances
that are near-optimally mapped on a given QPU. Our study mainly focuses on QA,
which are the most mature quantum technologies able to approximate large
optimization problems. However, our method can be extended to benchmark gate-
based QPU. The creation of the instance set is designed in a reverse fashion:
At first, a near-optimal mapping of a complete graph (source graph) on the
quantum chip (target graph) is found. The mapping function is then iteratively
altered to increase the size of the source graph while decreasing its density.
During this process, the target graph is not altered. The set of source graphs
is then used to generate instances of optimization problems with various
densities to benchmark the performance of QA. Our analysis provides valuable
insights into the characteristics of instances that could potentially lead to
a quantum advantage. Our finding is that despite near-optimal embeddings, the
first generation of D-Wave Advantage rapidly struggles to find good solutions
for instances with densities greater than 0.1.

The remainder of the paper is organized as follows: Section II introduces the
basic notions and notations on quantum annealing and minor-embedding. Section
III presents the method used to create the crafted instances and the technical
settings of the experiment. Section IV benchmarks the QPU against classical
solvers and discusses the results.

##  II Notation and definitions

Quantum annealers, such as D-Wave systems, constitute a noisy version of the
universal adiabatic quantum computer [23]. The evolution of these systems is
based on the interpolation of a Mixing Hamiltonian $H_{\mathrm{M}}$ whose
ground state is easy to prepare and a problem Hamiltonian $H_{\mathrm{P}}$
whose ground state encodes the solution to the problem. The problem
Hamiltonian can be fully specified with a source graph
$G_{\mathrm{s}}=(V_{\mathrm{s}},E_{\mathrm{s}})$:

$H_{\mathrm{P}}=\sum_{v\in V_{\mathrm{s}}}h_{v}\sigma_{v}^{z}+\sum_{(u,v)\in
E_{\mathrm{s}}}J_{uv}\sigma_{u}^{z}\sigma_{v}^{z}$ (1)

where $\sigma_{i}^{z}$ denotes the Pauli Z operator on variable $i$. The
ground state of $H_{\mathrm{p}}$ encodes the solution to the problem:

$\min C(\mathbf{s})=\sum_{v\in V_{\mathrm{s}}}h_{v}s_{v}+\sum_{(u,v)\in
E_{\mathrm{s}}}J_{uv}s_{u}s_{v}$ (2)

where $s_{i}\in\\{-1,+1\\}$. A Quadratic Unconstrained Binary Optimization
(QUBO) cost function with variables $x_{i}\in\\{0,1\\}$ is translated to an
Ising cost function 2 using a variable substitution $x_{i}=\frac{1+s_{i}}{2}$.
When a graph $G_{\mathrm{s}}$ cannot be straightly mapped to the qubits of the
quantum chip due to connectivity limitations, one has to find a mapping
function that maps the source graph $G_{\mathrm{s}}$ to a target graph
$G_{\mathrm{t}}=(V_{\mathrm{t}},E_{\mathrm{t}})$ where $G_{\mathrm{t}}$ is
defined as a subgraph of the quantum chip’s graph $G_{\mathrm{t\\_QPU}}$. This
problem is well-defined in the theory of graph minors [24] and is formally
stated as:

Find a function $\phi:V_{\mathrm{s}}\xrightarrow{}\mathcal{P}(V_{\mathrm{t}})$
such that :

  1. 1.

each vertex $v\in V_{\mathrm{s}}$ is mapped onto a connected subgraph
$\phi(v)$ of $G_{\mathrm{t}}$. $G_{\mathrm{t}}$ is a subgraph of
$G_{\mathrm{t\\_QPU}}$.

  2. 2.

each connected subgraph must be vertex disjoint
$\phi(v)\cap\phi(v^{\prime})=\emptyset$, with $v\neq v^{\prime}$.

  3. 3.

each edge $(u,v)\in E_{\mathrm{s}}$ is mapped onto at least one edge in
$E_{\mathrm{t}}$ : $\forall(u,v)\in E_{\mathrm{s}},\exists
u^{\prime}\in\phi(u),\exists
v^{\prime}\in\phi(v),~{}such~{}that~{}(u^{\prime},v^{\prime})\in
E_{\mathrm{t}}$.

Let $n_{\phi}$ be the number of nodes used in $G_{\mathrm{t}}$ to embed the
graph $G_{\mathrm{s}}$. A coupling strength is set to each edge of each
connected subgraph $\phi(v)$ to maintain ferromagnetic coupling between the
physical qubits, which act as a single logical variable. In this paper, we
refer to $v\in V_{s}$ as logical qubits or logical nodes (equivalent to
logical variables). We refer to $v\in V_{t}$ as physical qubits or physical
nodes. Each logical node is mapped to a set of physical nodes $\phi(v)$.

Two different types of algorithms are mainly used to find the function $\phi$.
The first set of algorithms takes as input both $G_{\mathrm{s}}$ and
$G_{\mathrm{t\\_QPU}}$. One implementation is the CMR heuristic [25], which
has the advantage of working with any graph $G_{\mathrm{t\\_QPU}}$. The second
set of algorithms take as input $G_{\mathrm{s}}$ and are designed for specific
$G_{\mathrm{t\\_QPU}}$. These algorithms usually take advantage of the
structure of $G_{\mathrm{t\\_QPU}}$ and produce mappings that are near-optimal
for complete graphs. For example, the CME heuristic [26] finds high-quality
mappings of complete graphs for D-Wave quantum chips, leveraging the regular
structure of the chip and considering inoperable qubits.

##  III Method

###  III-A Assessing the Quality of an Embedding

The minimum number of physical nodes required to embed a source graph
$G_{\mathrm{s}}$ can be lower bounded. Even if not reachable in practice, this
lower bound gives an idea of the quality of the embedding found and the
distance to an optimal embedding. Let $n_{\phi(v)^{*}}$ be a lower bound on
the optimal number of nodes required to embed a node $v\in V_{\mathrm{s}}$ on
the target graph. We assume $G_{\mathrm{t\\_QPU}}$ has a regular topology and
define $c_{\mathrm{phys}}$ as the number of edges per node, which is constant
for D-Wave topologies, _e.g._ , $c_{\mathrm{phys}}=15$ for Pegasus topology.
If $\phi(v)^{*}$ maps $v\in V_{\mathrm{s}}$ to a path (structure that
maximizes the potential connectivity of $\phi(v)^{*}$), the node $v$ requires
at least $n_{\phi(v)^{*}}$ nodes to be embedded on the target graph
$G_{\mathrm{t}}$:

$n_{\phi(v)^{*}}=\left.\begin{cases}1\text{ if }\mathit{deg}(v)\leq
c_{\mathrm{phys}}\\\ 2\text{ if
}c_{\mathrm{phys}}<\mathit{deg}(v)\leq(2c_{\mathrm{phys}}-2)\\\
\left\lceil\frac{\mathit{deg}(v)-(2c_{\mathrm{phys}}-2)}{c_{\mathrm{phys}}-2}\right\rceil+2\text{
otherwise }\end{cases}\right\\}.$ (3)

The number of nodes used in the target graph to embed $G_{\mathrm{s}}$ is then
$n_{\phi}\geq\sum_{v\in V_{\mathrm{s}}}n_{\phi(v)^{*}}$. This bound is then
used to define an overhead ratio $r_{\mathrm{o}}$ that computes the overhead
of physical qubits used by the mapping function:

$r_{\mathrm{o}}=\frac{n_{\phi}}{\sum_{v\in V_{s}}n_{\phi(v)^{*}}}.$ (4)

The closer the ratio $r_{\mathrm{o}}$ is to 1, the better. We use this metric
to show that our instance generation process creates near-optimally embedded
instances.

###  III-B Near-optimal Mapping Generation

Fig. 1 shows the workflow creating graph instances near-optimally embedded on
the D-Wave quantum chip. The first step consists of finding the mapping
function of a complete graph $G_{\mathrm{s}}$ to a target graph
$G_{\mathrm{t}}$ that is a subgraph of $G_{\mathrm{t\\_QPU}}$. The CME
algorithm finds such mapping (see Fig. 1a). Each logical node $v\in
V_{\mathrm{s}}$ is mapped to a connected subgraph $\phi_{\mathrm{CME}}(v)$ on
the D-Wave chip, which topology is represented by the graph
$G_{\mathrm{t\\_QPU}}$. One special feature is that the CME algorithm
generates subgraphs $\phi_{\mathrm{CME}}(v)$ that are paths to maximize the
potential connectivity of each logical qubit $\phi_{\mathrm{CME}}(v)$. It is
then possible to increase the size of the logical graph $G_{\mathrm{s}}$ by
selecting and splitting a random path $\phi_{\mathrm{CME}}(v)$ in two equal
parts (see Fig. 1b and c). The new source graph
$G_{\mathrm{s}}^{\prime}=(V_{\mathrm{s}}^{\prime},E_{\mathrm{s}}^{\prime})$
has a lower density than $G_{\mathrm{s}}$, but an additional node owing to the
split. This last step is repeated to create arbitrarily dense source graphs
$G_{\mathrm{s}}^{\prime}$, preserving the near-optimal embedding clause in the
new mapping function $\phi_{\mathrm{nopt}}(v)$. As the number of iterations
increases, the Graph Edit Distance (GED) between $G_{\mathrm{s}}^{\prime}$ and
$G_{\mathrm{t}}$ shrinks, meaning that the topology of
$G_{\mathrm{s}}^{\prime}$ becomes closer to the topology of $G_{\mathrm{t}}$.
It is important to notice that the set of nodes $V_{\mathrm{t}}$ of
$G_{\mathrm{t}}$ are defined by the CME method at the beginning and remain
fixed during the steps used to generate the different instances
$G_{\mathrm{s}}^{\prime}$. Hence, all the instances $G_{\mathrm{s}}^{\prime}$
use the same physical graph $G_{\mathrm{t}}$. Only the mapping function
$\phi_{\mathrm{nopt}}$ varies between each graph $G_{\mathrm{s}}^{\prime}$.
The generated graphs $G_{\mathrm{s}}^{\prime}$ that are very sparse are
considered favorable to the quantum computer as they have many logical
variables, each mapped to very few or only one physical qubit. On the
contrary, graphs $G_{\mathrm{s}}^{\prime}$ with high densities have each
vertex $v$ mapped to long chains of physical qubits $\phi_{\mathrm{nopt}}(v)$.
It suggests that this type of instance is less favorable to the quantum
annealer and will be more easily solved by classical computers.

Figure 1: Instance generation with iterative chain split. a) The algorithm
starts with a complete graph embedding using CME method on the graph
$G_{\mathrm{t\\_QPU}}$. Each node in the source graph $G_{\mathrm{s}}$ is
mapped to a chain in the target graph $G_{\mathrm{t}}$. b) A random logical
node $v\in V_{\mathrm{s}}$ is selected. c) The corresponding chain
$\phi_{\mathrm{CME}}(v)$ is split in two parts to create a new logical node in
the source graph. It changes the mapping function to $\phi_{\mathrm{nopt}}$.
Steps b. and c. are repeated on $G_{\mathrm{s}}^{\prime}$ until the desired
density in the logical graph is reached.

We evaluate the efficiency of our method using a complete graph of size
$n=100$. The CME method finds the mapping function of the complete graph
$G_{\mathrm{s}}$ to D-Wave Advantage6.4, using 982 physical qubits. The method
described in the previous section is used to generate source graphs
$G_{\mathrm{s}}^{\prime}$ of diverse densities $d\in\\{0.9,0.8,...,0.1\\}$
with their associated near-optimal mapping $\phi_{\mathrm{nopt}}$. For each
source graph $G_{\mathrm{s}}^{\prime}$, the CMR method is run 100 times with a
time limit set to 5 minutes to find a second mapping
$\phi_{\mathrm{100\\_CMR}}$ that uses the least number of nodes in
$G_{\mathrm{t\\_QPU}}$. Fig. 2a compares the values of $r_{\mathrm{o}}$ for
mappings: $\phi_{\mathrm{nopt}}$ and $\phi_{\mathrm{100\\_CMR}}$ and clearly
shows that $\phi_{\mathrm{nopt}}$ mappings use a number of nodes close to the
optimal bound defined in 3. The mappings $\phi_{\mathrm{100\\_CMR}}$ use more
nodes that $\phi_{\mathrm{nopt}}$ for every instances. The worst case behavior
of $\phi_{\mathrm{100\\_CMR}}$ is on sparse instances, whereas
$\phi_{\mathrm{nopt}}$ produces embeddings that can be considered near-
optimal. We then use the graphs $G_{\mathrm{s}}^{\prime}$ to create max-cut
instances and compute the ratio of cut sizes obtained using each mapping
$\phi_{\mathrm{nopt}}$ and $\phi_{\mathrm{100\\_CMR}}$ with D-Wave
Advantage6.4. Fig. 2b shows that the best-cut size is increased when using the
mapping $\phi_{\mathrm{nopt}}$ instead of the mapping
$\phi_{\mathrm{100\\_CMR}}$. It confirms that the mappings found by our
method, in addition to using fewer qubits, increase the performance of the
D-Wave Advantage6.4.

###  III-C Instance Set

In the rest of the paper, we generate instances from a complete graph of size
$n=174$. It is the largest clique that can be embedded with the CME method on
Advantage6.4, using $2918$ physical qubits, representing approximately $52\%$
of the physical working qubits of the chip. Using the method described above,
we generate graphs of various densities $G_{\mathrm{s}}^{\prime}$ with their
associated near-optimal transformation $\phi_{\mathrm{nopt}}$ on Advantage6.4
(see Table II for instances description). The source graphs
$G_{\mathrm{s}}^{\prime}$ are then used to create instances of three different
optimization problems with 30 instances for each density. The two first
problems are unweighted and weighted max-cut problems defined as:

Figure 2: Performance comparison of mapping functions found by our method
$\phi_{\mathrm{nopt}}$ and the best among 100 tries of CMR method
$\phi_{\mathrm{100\\_CMR}}$. 30 instances are generated for each density. a)
Comparison using the overhead ratio (see 4). Error bars show the standard
deviations b) Ratio of the best cut size obtained for each mapping for each
instance.

$\min C(\mathbf{s})=\sum_{(u,v)\in E_{\mathrm{s}}^{\prime}}J_{uv}s_{u}s_{v}$
(5)

where $s_{u},s_{v}\in\\{-1,+1\\}$. $J_{uv}=1$ for unweighted max-cut and
$J_{uv}\in\\{-\frac{128}{128},...,-\frac{1}{128},\frac{1}{128},...,\frac{128}{128}\\}$
for weighted max-cut. Weights are used to challenge the precision and long-
range coherence of the QPU. The precision is the same as in [4]. The third
problem is a weighted Maximum Independent Set (MIS) defined by the cost
function:

$\min C(\mathbf{x})=-\sum_{v\in
V_{\mathrm{s}}^{\prime}}\omega_{v}x_{v}+\sum_{(u,v)\in
E_{\mathrm{s}}^{\prime}}\omega_{uv}x_{u}x_{v}$ (6)

where $x_{u},x_{v}\in\\{0,1\\}$,
$\omega_{v}\in\\{\frac{1}{128},...,\frac{128}{128}\\}$ and $\omega_{uv}=2$.
$\omega_{uv}$ acts as the penalty term used to enforce the hard constraints
for D-Wave and Tabu search. The constraints are encoded with inequalities on
the solver Gurobi to reduce the search space.

###  III-D Solver Settings

####  III-D1 Advantage6.4 settings

The performance of the quantum annealer is evaluated using an access time
limit of 1 second. Each instance is mapped on the QPU using the mapping
$\phi_{\mathrm{nopt}}$. The repartition of the processing time considering
different annealing times $\\{1,10,100,1000\\}\mu s$ with respectively
$\\{5000,4780,3310,820\\}$ shots is illustrated in Fig. 3. It shows the
existing trade-offs between the number of shots and the duration of the
annealing time per shot. Fig. 4 shows the performance of D-Wave considering
different values of annealing time for each problem. The average optimal
annealing time seems to be problem-dependent. Considering these results, we
choose $1\mu s$ and $5000$ shots for unweighted max-cut, $1000\mu s$ and $820$
shots for weighted max-cut and $1\mu s$ and $5000$ shots for weighted MIS
problems. The uniform_torque_compensation method implemented on D-Wave sets
the chain strength used to maintain ferromagnetic couplings between the
physical qubits $\phi_{\mathrm{nopt}}(v)$. The weight associated with a single
qubit $h_{v}$ is uniformly spread over the physical qubits
$\phi_{\mathrm{nopt}}(v)$. The same method is used for edges. A majority vote
is used to unembed the problem. We do not use annealing offset or spin
reversal methods. MIS instances are post-processed to avoid constraint
violations. We iteratively and randomly select edges that violate the
constraint and only keep the node with the higher weight in the independent
set. The benchmark is done using the best solution found among all shots for
each instance.

####  III-D2 Tabu Search settings

The Tabu search is based on the method introduced in [27]. The Tabu search is
run directly on the instance associated with the source graph
$G_{\mathrm{s}}^{\prime}$. The processing time limit is set to $1s$ and the
tenure length is set to $|V_{\mathrm{s}}^{\prime}|/4$. The Tabu search uses a
single core of processor AMD EPYC 7702P. MIS instances are post-processed with
the same method used for D-Wave.

TABLE II: Properties of $G_{\mathrm{s}}^{\prime}$ graphs used to build optimization problem instances. These graphs are mapped on 2918 qubits. Density | Avg $|V_{\mathrm{s}}^{\prime}|$ | $r_{\mathrm{o}}$ | Density | Avg $|V_{\mathrm{s}}^{\prime}|$ | $r_{\mathrm{o}}$
---|---|---|---|---|---
0.02 | 1318 | 1.06 | 0.1 | 565 | 1.14
0.03 | 1062 | 1.08 | 0.2 | 395 | 1.16
0.04 | 912 | 1.10 | 0.3 | 321 | 1.17
0.05 | 810 | 1.11 | 0.4 | 277 | 1.18
0.06 | 737 | 1.12 | 0.5 | 248 | 1.19
0.07 | 680 | 1.12 | 0.6 | 226 | 1.19
0.08 | 635 | 1.13 | 0.7 | 209 | 1.19
0.09 | 597 | 1.13 | 0.8 | 195 | 1.20
|  |  | 0.9 | 184 | 1.20

Figure 3: QPU access time repartition for different annealing times per shot
with a total running time limit of $1s$. Each color represents the fraction of
time used for programming the QPU, annealing the $n$ shots, reading the
results of the $n$ shots, and delays between the $n$ shots.

####  III-D3 Gurobi settings

Gurobi [28] is run with a time limit of $60s$. It is used to establish a
reference solution with which the Tabu search and D-Wave are compared. The
solver is set with default parameters. The solver is run on the instance
associated with the source graph $G_{\mathrm{s}}^{\prime}$. For the MIS
problem, the constraints are directly encoded with inequalities. Gurobi is
parallelized on 20 cores of a processor AMD EPYC 7702P.

####  III-D4 Random solver settings

A random solver is used on MIS instances. It iteratively and randomly selects
a node in the set $V_{\mathrm{s}}^{\prime}$ and removes from
$V_{\mathrm{s}}^{\prime}$ all the neighbors of this node. This method
generates the same number of solutions generated by D-Wave, i.e., 5000
solutions. The best MIS is then selected among these shots. We did not use
this method for the max-cut problem as random solution cuts are far from being
competitive with D-Wave and Tabu search solutions.

Figure 4: Annealing time scan with a quantum processing time limit of $1s$ for
three types of problems with various densities. a) Unweighted maxcut instances
b) Weighted maxcut instances c) Weighted MIS.

##  IV Results and Discussion

We compare the performance of D-Wave Advantage6.4 against classical solvers to
approximate solutions of optimization problems introduced in the previous
section. Table II shows the average number of logical variables for each
density processed in Fig. 5. We recall that the density of the graph is
inversely proportional to the number of variables in the graph
$G_{s}^{\prime}$. Hence, the size of instances increases when the density
decreases. The results of D-Wave and the Tabu search are expressed as a ratio
of a reference solution computed with Gurobi with a runtime limit set to
$60s$. In general, D-Wave performs well with respect to Tabu search on sparse
instances under 0.1 density. Under this density, the logical problem has more
than 565 variables, which limits the performance of the Tabu search due to the
large size of the solution space. We observed that the Tabu search did not
have sufficient processing time for large instances, leading to almost random
solutions to these instances. It explains why the box plots are stretched for
Tabu search for small densities. The performance transition between the
quantum and classical heuristics occurs at density $d=0.06$ for unweighted
max-cut, around $d=0.1$ for weighted max-cut instances and at $d=0.05$ for MIS
instances. On very sparse instances ($d<0.03$), the Advantage6.4 is able to
outperform Gurobi on few instances of unweighted and weighted max-cut, showing
that the QPU may be useful against classical methods for very sparse instances
having more than 1000 variables. For denser instances, lots of extra qubits
are used to encode the problem on the QPU. These instances become easier for
classical solvers as the search space reduces exponentially with the density
but stays approximately the same for the QPU as the number of physical qubits
remains the same due to the embedding. In addition, long chains of physical
qubits require strong ferromagnetic couplings to maintain spin orders in
chains $\phi_{\mathrm{nopt}}(v)$. It adds extra energy to the system and
impacts the couplers’ precision due to the rescaling of Hamiltonian weights.
Both Tabu search and D-Wave QPU do not seem to be adapted to approximate
solutions to the MIS problem. These methods are unable to compete with the
random solver used to generate valid maximal independent sets. D-Wave only
starts to be competitive on instances with densities lower than 0.02. This
behavior is reminiscent of the fact that both D-Wave and Tabu search
performances are limited by the constraints encoding using penalty methods and
cannot produce valid solutions in a dense regime.

Figure 5: Benchmark of D-Wave Advantage6.4 and Tabu search (1s time limit).
The results are expressed as ratio of the best result found by Gurobi (60s
runtime limit). The size of each instance is provided in Table II a)
Unweighted max-cut b) Weighted max-cut c) MIS problem with a random solver.

Several points have to be clarified and considered carefully to avoid
misleading conclusions. The aim of this benchmark is not to study the scaling
of the QPU but rather to determine when the quantum computer performs best at
solving instances that use the same physical requirements (in our case, the
number of physical qubits). The settings of this experiment only involved
tuning the annealing time, which has a huge impact when working on time-
limited tasks. We did not want to include other complex parameter optimization
such as the advanced tuning of the chain strength [29, 30, 31], annealing
schedules [32, 33] or advanced calibration methods called shimming [34]. The
main idea of our study is to provide an objective view of D-Wave’s average
performance without employing complex preprocessing methods. Nevertheless, the
study of the trade-offs brought by these methods represents a relevant
perspective for future work. The comparison with classical heuristics can also
be enhanced by using algorithms that excel on sparse instances of the maxcut
problem [35].

The minor-embedding task is a step that is usually excluded from quantum
benchmark studies, which almost only consider the pure annealing time without
considering delays, readouts and programming time. However, a huge amount of
classical preprocessing time can be spent to optimize the mapping function,
making comparing existing studies difficult. Instead, our approach permits the
evaluation of the QPU on near-optimally mapped instances, giving insights into
the properties that leverage the power of QA.

To conclude the discussion, this benchmark methodology is not fair and
representative of the performance of D-Wave systems’ ability to solve real-
world instances. It rather acts as a performance evaluation of D-Wave systems’
ability to solve large optimization instances based on an ideal mapping of the
instance on the QPU, using QA default settings. This methodology could be used
to detect optimization problems that could benefit from a quantum advantage.

##  V Conclusion

This article introduced a new method to create source graphs of various
densities near-optimally mapped on a QA. This set of graphs was then used to
generate large optimization problem instances to benchmark D-Wave QPU against
classical solvers. Considering an equal processing time limit, the QPU was
competitive with Tabu search for instances with densities inferior to $0.1$ on
the Weighted max-cut problem. It even outperformed Gurobi solver on some
really sparse unweighted max-cut instances, whereas Gurobi had $60$ times the
processing time of D-Wave QPU. Conversely, the QPU seemed less efficient at
approximating similar-sized constrained problems. It is important to remember
that this set of crafted instances is very favorable and specifically designed
for D-Wave QA. Our results can be seen as the benchmark of the default
behavior of D-Wave QPU on large instances of ideal shape. On real instances,
the mapping function found by embedding methods such as CMR will not be as
efficient as the one designed in this paper.

A relevant perspective is to study other optimization problems using this
methodology to build the instance set. It could also be used to evaluate the
benefits of advanced preprocessing methods of quantum annealers. Finally, this
methodology can also be adapted to benchmark variational quantum algorithms
working on gate-based quantum computers, considering the generation of
instances adapted to the topology of the quantum chip with a planted optimal
swapping network.

## References

  * [1] F. Arute, K. Arya, R. Babbush, et al., “Quantum supremacy using a programmable superconducting processor,” Nature, vol. 574, no. 7779, p. 505–510, 2019.
  * [2] H.-S. Zhong, H. Wang, Y.-H. Deng, et al., “Quantum computational advantage using photons,” Science, vol. 370, no. 6523, pp. 1460–1463, 2020\.
  * [3] L. S. Madsen, F. Laudenbach, M. F. Askarani, et al., “Quantum computational advantage with a programmable photonic processor,” Nature, vol. 606, no. 7912, p. 75–81, 2022.
  * [4] A. D. King, A. Nocera, M. M. Rams, et al., “Computational supremacy in quantum simulation,” arXiv preprint arXiv:2403.00910, 2024.
  * [5] S. Martiel, T. Ayral, and C. Allouche, “Benchmarking quantum coprocessors in an application-centric, hardware-agnostic, and scalable way,” IEEE Transactions on Quantum Engineering, vol. 2, pp. 1–11, 2021.
  * [6] D. Mills, S. Sivarajah, T. L. Scholten, and R. Duncan, “Application-motivated, holistic benchmarking of a full quantum computing stack,” Quantum, vol. 5, p. 415, 2021.
  * [7] J. R. Finžgar, P. Ross, L. Hölscher, et al., “Quark: A framework for quantum computing application benchmarking,” in 2022 IEEE international conference on quantum computing and engineering (QCE), pp. 226–237, IEEE, 2022.
  * [8] T. Lubinski, S. Johri, P. Varosy, et al., “Application-oriented performance benchmarks for quantum computing,” IEEE Transactions on Quantum Engineering, 2023.
  * [9] F. Barbaresco, L. Rioux, C. Labreuche, et al., “Bacq–application-oriented benchmarks for quantum computing,” arXiv preprint arXiv:2403.12205, 2024.
  * [10] T. Tomesh, P. Gokhale, V. Omole, et al., “Supermarq: A scalable quantum benchmark suite,” in 2022 IEEE International Symposium on High-Performance Computer Architecture (HPCA), pp. 587–603, IEEE, 2022.
  * [11] V. Gilbert, S. Louise, and R. Sirdey, TAQOS: A Benchmark Protocol for Quantum Optimization Systems, p. 168–176.  Springer Nature Switzerland, 2023.
  * [12] N. P. Sawaya, D. Marti-Dafcik, Y. Ho, et al., “Hamlib: A library of hamiltonians for benchmarking quantum algorithms and hardware,” in 2023 IEEE International Conference on Quantum Computing and Engineering (QCE), vol. 2, pp. 389–390, IEEE, 2023.
  * [13] N. Quetschlich, L. Burgholzer, and R. Wille, “Mqt bench: Benchmarking software and design automation tools for quantum computing,” Quantum, vol. 7, p. 1062, 2023.
  * [14] I. Dunning, S. Gupta, and J. Silberholz, “What works best when? a systematic evaluation of heuristics for max-cut and qubo,” INFORMS Journal on Computing, vol. 30, no. 3, p. 608–624, 2018.
  * [15] F. Furini, E. Traversi, P. Belotti, et al., “Qplib: a library of quadratic programming instances,” Mathematical Programming Computation, vol. 11, no. 2, p. 237–265, 2018.
  * [16] A. Abbas, A. Ambainis, B. Augustino, et al., “Quantum optimization: Potential, challenges, and the path forward,” arXiv preprint arXiv:2312.02279, 2023.
  * [17] “D-wave system. solver properties and parameters.”  https://docs.dwavesys.com/docs/latest/doc_solver_ref.html [Accessed 20-04-2024].
  * [18] V. Choi, “Minor-embedding in adiabatic quantum computation: I. the parameter setting problem,” Quantum Information Processing, vol. 7, pp. 193–209, 2008.
  * [19] E. Pelofske, A. Bärtschi, and S. Eidenbenz, “Quantum annealing vs. qaoa: 127 qubit higher-order ising problems on nisq computers,” in International Conference on High Performance Computing, pp. 240–258, Springer, 2023.
  * [20] Y. Pang, C. Coffrin, A. Y. Lokhov, and M. Vuffray, “The potential of quantum annealing for rapid solution structure identification,” Constraints, vol. 26, no. 1, pp. 1–25, 2021.
  * [21] B. Tasseff, T. Albash, Z. Morrell, et al., “On the emerging potential of quantum annealing hardware for combinatorial optimization,” arXiv preprint arXiv:2210.04291, 2022.
  * [22] T. Lubinski, C. Coffrin, C. McGeoch, et al., “Optimization applications as quantum performance benchmarks,” arXiv preprint arXiv:2302.02278, 2023\.
  * [23] T. Albash and D. A. Lidar, “Adiabatic quantum computation,” Reviews of Modern Physics, vol. 90, no. 1, p. 015002, 2018.
  * [24] N. Robertson and P. Seymour, “Graph minors .xiii. the disjoint paths problem,” Journal of Combinatorial Theory, Series B, vol. 63, no. 1, p. 65–110, 1995.
  * [25] J. Cai, W. G. Macready, and A. Roy, “A practical heuristic for finding graph minors,” arXiv preprint arXiv:1406.2741, 2014.
  * [26] T. Boothby, A. D. King, and A. Roy, “Fast clique minor generation in chimera qubit connectivity graphs,” Quantum Information Processing, vol. 15, pp. 495–508, 2016.
  * [27] G. Palubeckis, “Multistart tabu search strategies for the unconstrained binary quadratic optimization problem,” Annals of Operations Research, vol. 131, no. 1–4, p. 259–282, 2004.
  * [28] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2023.
  * [29] V. Gilbert and S. Louise, “Quantum annealers chain strengths: A simple heuristic to set them all,” arXiv preprint arXiv:2404.05443, 2024.
  * [30] D. Willsch, M. Willsch, C. D. Gonzalez Calaza, et al., “Benchmarking advantage and d-wave 2000q quantum annealers with exact cover problems,” Quantum Information Processing, vol. 21, no. 4, p. 141, 2022.
  * [31] T. V. Le, M. V. Nguyen, T. N. Nguyen, et al., “Benchmarking chain strength: An optimal approach for quantum annealing,” in 2023 IEEE International Conference on Quantum Computing and Engineering (QCE), IEEE, 2023\.
  * [32] J. I. Adame and P. L. McMahon, “Inhomogeneous driving in quantum annealers can result in orders-of-magnitude improvements in performance,” Quantum Science and Technology, vol. 5, no. 3, p. 035011, 2020.
  * [33] M. Khezri, X. Dai, R. Yang, et al., “Customized quantum annealing schedules,” Physical Review Applied, vol. 17, no. 4, 2022.
  * [34] K. Chern, K. Boothby, J. Raymond, et al., “Tutorial: Calibration refinement in quantum annealing,” arXiv preprint arXiv:2304.10352, 2023\.
  * [35] D. Rehfeldt, T. Koch, and Y. Shinano, “Faster exact solution of sparse maxcut and qubo problems,” Mathematical Programming Computation, vol. 15, no. 3, p. 445–470, 2023.