EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-PH-EP-2015-098 LHCb-PAPER-2015-015 30 June 2015

Quantum numbers of the $X(3872)$ state and orbital angular momentum in its
$\rho^{0}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ decay

The LHCb collaboration†††Authors are listed at the end of this paper.

Angular correlations in $B^{+}\rightarrow X(3872)K^{+}$ decays, with
$X(3872)\rightarrow\rho^{0}J/\psi$, $\rho^{0}\rightarrow\pi^{+}\pi^{-}$ and
$J/\psi\rightarrow\mu^{+}\mu^{-}$, are used to measure orbital angular
momentum contributions and to determine the $J^{PC}$ value of the $X(3872)$
meson. The data correspond to an integrated luminosity of 3.0 fb-1 of proton-
proton collisions collected with the LHCb detector. This determination, for
the first time performed without assuming a value for the orbital angular
momentum, confirms the quantum numbers to be $J^{PC}=1^{++}$. The $X(3872)$ is
found to decay predominantly through S wave and an upper limit of $4\%$ at
$95\%$ C.L. is set on the D–wave contribution.

Published in Phys. Rev. D92, 011102 (2015).

© CERN on behalf of the LHCb collaboration, license CC-BY-4.0.

The $X(3872)$ state was discovered in $B^{+,0}\rightarrow X(3872)K^{+,0}$,
$X(3872)\rightarrow\pi^{+}\pi^{-}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}$, ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}\rightarrow\ell^{+}\ell^{-}$ decays by the Belle experiment [1] and
subsequently confirmed by other experiments [2, 3, 4].111 The inclusion of
charge-conjugate states is implied in this article. Its production was also
studied at the LHC [5, 6]. However, the nature of this state remains unclear.
The $X(3872)$ state is narrow, has a mass very close to the $D^{0}\kern
1.99997pt\overline{\kern-1.99997ptD}{}^{*0}$ threshold and decays to
$\rho^{0}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ and
$\omega{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ final states with
comparable branching fractions [7], thus violating isospin symmetry. This
suggests that the $X(3872)$ particle may not be a simple $c\bar{c}$ state, and
exotic states such as $D^{0}\kern 1.99997pt\overline{\kern-1.99997ptD}{}^{*0}$
molecules [8], tetraquarks [9] or mixtures of states [10] have been proposed
to explain its composition. The $X(3872)$ quantum numbers, such as total
angular momentum $J$, parity $P$, and charge-conjugation $C$, impose
constraints on the theoretical models of this state. The orbital angular
momentum $L$ in the $X(3872)$ decay may also provide information on its
internal structure.

Observations of the
$X(3872)\rightarrow\gamma{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ and
$X(3872)\rightarrow\gamma\psi(2S)$ decays [11, 12, 13] imply positive $C$,
which requires the total angular momentum of the dipion system ($J_{\pi\pi}$)
in $X(3872)\rightarrow\pi^{+}\pi^{-}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}$ decays to be odd. Dipion mass, $M(\pi^{+}\pi^{-})$, is limited by the
available phase space to be less than $775$ MeV, and so $J_{\pi\pi}\geq 3$ can
be ruled out since there are no known or predicted mesons with such high spins
at such low masses.222 We use mass and momentum units in which $c=1$. In fact,
the distribution of $M(\pi^{+}\pi^{-})$ is consistent with
$X(3872)\rightarrow\rho^{0}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$
decays [14, 15, 6], in line with $J_{\pi\pi}=1$, the only plausible value.

The choices for $J^{PC}$ were narrowed down to two possibilities, $1^{++}$ or
$2^{-+}$, by the CDF collaboration, via an analysis of the angular
correlations in inclusively reconstructed
$X(3872)\rightarrow\pi^{+}\pi^{-}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$
and ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}\rightarrow\mu^{+}\mu^{-}$
decays, dominated by prompt production in $p\bar{p}$ collisions [16]. Using
1.0 fb-1 of $pp$ collision data collected by LHCb, $J^{PC}=2^{-+}$ was ruled
out in favor of the $1^{++}$ assignment, using the angular correlations in the
same decay chain, with the $X(3872)$ state produced in $B^{+}\rightarrow
X(3872)K^{+}$ decays [17]. Both angular analyses assumed that the lowest
orbital angular momentum between the $X(3872)$ decay products ($L_{\rm min}$)
dominated the matrix element. Significant contributions from $L_{\rm min}+2$
amplitudes could invalidate the $1^{++}$ assignment. Since the phase-space
limit on $M(\pi^{+}\pi^{-})$ is close to the $\rho^{0}$ pole ($775.3\pm 0.3$
MeV [7]), the energy release in the $X(3872)$ decay, $Q\equiv
M({J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}\pi^{+}\pi^{-})-M({J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu})-M(\pi^{+}\pi^{-})$, is a small fraction of the $X(3872)$ mass, making
the orbital angular momentum barrier effective.333 Dimuon candidates are
constrained to the known ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ mass
[7]. However, an exotic component in $X(3872)$ could induce contributions from
higher orbital angular momentum for models in which the size of the $X(3872)$
state is substantially larger than the compact sizes of the charmonium states.
Therefore, it is important to probe the $X(3872)$ spin-parity without any
assumptions about $L$. A determination of the magnitude of contributions from
$L_{\rm min}+2$ amplitudes for the correct $J^{PC}$ is also of interest, since
a substantial value would suggest an anomalously large size of the $X(3872)$
state. In this article, we extend our previous analysis [17] of five-
dimensional angular correlations in $B^{+}\rightarrow X(3872)K^{+}$,
$X(3872)\rightarrow\rho^{0}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$,
$\rho^{0}\rightarrow\pi^{+}\pi^{-}$, ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}\rightarrow\mu^{+}\mu^{-}$ decays to accomplish these goals. The
integrated luminosity of the data sample has been tripled by adding $8$ TeV
$pp$ collision data collected in 2012.

The LHCb detector is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$, described in detail in Ref. [18, 19]. The
$X(3872)$ candidate selection, which is based on reconstructing
$B^{+}\rightarrow({J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}\rightarrow\mu^{+}\mu^{-})\pi^{+}\pi^{-}K^{+}$ candidates using particle
identification information, transverse momentum ($p_{\rm T}$) thresholds and
requiring separation of tracks and the $B^{+}$ vertex from the primary $pp$
interaction vertex, is improved relative to that of Ref. [17]. The signal
efficiency is increased by lowering requirements on $p_{\rm T}$ for muons from
0.90 to 0.55 GeV and for hadrons from 0.25 to 0.20 GeV. The background is
further suppressed without significant loss of signal by requiring $Q<250$
MeV. The $X(3872)$ mass resolution ($\sigma_{\Delta M}$) is improved from
about $5.5$ MeV to $2.8$ MeV by constraining the $B^{+}$ candidate to its
known mass and requiring its momentum to point to a $pp$ collision vertex in
the kinematic fit of its decay. The distribution of $\Delta M\equiv
M(\pi^{+}\pi^{-}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu})-M({J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu})$ is shown in Fig. 1. A
Crystal Ball function [20] with symmetric tails is used to model the signal
shape, while the background is assumed to be linear. An unbinned maximum
likelihood fit yields $1011\pm 38$ $B^{+}\rightarrow X(3872)K^{+}$ decays and
$1468\pm 44$ background entries in the $725<\Delta M<825$ MeV range used in
the angular analysis. The signal purity is $80\%$ within $2.5\,\sigma_{\Delta
M}$ from the signal peak. From studying the $K^{+}\pi^{+}\pi^{-}$ mass
distribution, the dominant source of background is found to be
$B^{+}\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K_{1}(1270)^{+}$,
$K_{1}(1270)^{+}\rightarrow K^{+}\pi^{+}\pi^{-}$ decays.

Figure 1: Distribution of $\Delta M$ for
$B^{+}\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}K^{+}\pi^{+}\pi^{-}$ candidates. The fit of the $X(3872)$ signal is
displayed. The solid (blue), dashed (red), and dotted (green) lines represent
the total fit, signal component, and background component, respectively.

Angular correlations in the $B^{+}$ decay chain are analyzed using an unbinned
maximum-likelihood fit to determine the $X(3872)$ quantum numbers and orbital
angular momentum in its decay. The probability density function (${\cal P}$)
for each $J^{PC}$ hypothesis, $J_{X}$, is defined in the five-dimensional
angular space $\Omega\equiv$
$(\cos\theta_{X},\cos\theta_{\rho},\Delta\phi_{X,\rho},\cos\theta_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}},\Delta\phi_{X,{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}})$, where
$\theta_{X}$, $\theta_{\rho}$ and
$\theta_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}$ are the helicity
angles [21, 22, 23] in $X(3872)$, $\rho^{0}$ and
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ decays, respectively, and
$\Delta\phi_{X,\rho}$, $\Delta\phi_{X,{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}$ are the angles between the decay planes of the $X(3872)$ particle and
of its decay products. The quantity ${\cal P}$ is the normalized product of
the expected decay matrix element (${\cal M}$) squared and of the
reconstruction efficiency ($\epsilon$), ${\cal P}(\Omega|J_{X})=|{\cal
M}(\Omega|J_{X})|^{2}\,\epsilon(\Omega)/I(J_{X})$, where $I(J_{X})=\int|{\cal
M}(\Omega|J_{X})|^{2}\,\epsilon(\Omega){\it d}\Omega$. The efficiency is
averaged over the $\pi^{+}\pi^{-}$ mass using a simulation [24, 25, 26, 27,
*Agostinelli:2002hh, 29] of the
$X(3872)\rightarrow\rho^{0}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$,
$\rho^{0}\rightarrow\pi^{+}\pi^{-}$ decay. The lineshape of the $\rho^{0}$
resonance can change slightly depending on the $X(3872)$ spin hypothesis. The
effect on $\epsilon(\Omega)$ is very small and is neglected. The angular
correlations are obtained using the helicity formalism [16],

$\displaystyle|\,{\cal M}(\Omega|J_{X})\,|^{2}$ $\displaystyle=$
$\displaystyle\sum_{\Delta\lambda_{\mu}=-1,+1}\left|\phantom{A_{\lambda_{\psi}}}\hskip-11.38092pt\right.\sum_{\lambda_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}},\lambda_{\rho}=-1,0,+1}A_{\lambda_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}},\lambda_{\rho}}\quad
D^{J_{X}}_{0\,,\,\lambda_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}-\lambda_{\rho}}(0,\theta_{X},0)^{*}$ $\displaystyle\hskip
68.28644ptD^{1}_{\lambda_{\rho}\,,\,0}(\Delta\phi_{X,\rho},\theta_{\rho},0)^{*}\quad
D^{1}_{\lambda_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}\,,\,\Delta\lambda_{\mu}}(\Delta\phi_{X,{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}},\theta_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}},0)^{*}\left.\phantom{A_{\lambda_{\psi}}}\hskip-11.38092pt\right|^{2},$

where $\lambda$ are particle helicities,
$\Delta\lambda_{\mu}=\lambda_{\mu^{+}}-\lambda_{\mu^{-}}$, and
$D^{J}_{\lambda_{1}\,,\,\lambda_{2}}$ are Wigner functions [21, 22, 23]. The
helicity couplings, $A_{\lambda_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}},\lambda_{\rho}}$, are expressed in terms of the $LS$ couplings,
$B_{LS}$, with the help of Clebsch-Gordan coefficients, where $L$ is the
orbital angular momentum between the $\rho^{0}$ and the
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ mesons, and $S$ is the sum of
their spins,

$\displaystyle A_{\lambda_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}\,,\,\lambda_{\rho}}=\sum_{L}\sum_{S}B_{LS}\left(\begin{array}[]{cc|c}J_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}&J_{\rho}&S\\\ \lambda_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}&-\lambda_{\rho}&\lambda_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}-\lambda_{\rho}\end{array}\right)\left(\begin{array}[]{cc|c}L&S&J_{X}\\\
0&\lambda_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}-\lambda_{\rho}&\lambda_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}-\lambda_{\rho}\end{array}\right).$ (6) (7)

Possible values of $L$ are constrained by parity conservation,
$P_{X}=P_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}\,P_{\rho}\,(-1)^{L}=(-1)^{L}$. In the previous analyses [16, 14, 17],
only the minimal value of the angular momentum, $L_{\rm min}$, was allowed.
Thus, for the preferred $J^{PC}=1^{++}$ hypothesis, the D wave was neglected
allowing only S–wave decays. In this work all $L$ values are allowed in Eq.
(7). The corresponding $B_{LS}$ amplitudes are listed in Table 1. Values of
$J_{X}$ up to four are analyzed. Since the orbital angular momentum in the
$B^{+}$ decay equals $J_{X}$, high values are suppressed by the angular
momentum barrier. In fact, the highest observed spin of any resonance produced
in $B$ decays is three [30, 31]. Since ${\cal P}$ is insensitive to the
overall normalization of the $B_{LS}$ couplings and to the phase of the matrix
element, the $B_{LS}$ amplitude with the lowest $L$ and $S$ is set to the
arbitrary reference value $(1,0)$. The set of other possible complex $B_{LS}$
amplitudes, which are free parameters in the fit, is denoted as $\alpha$.

Table 1: Parity-allowed $LS$ couplings in the $X(3872)\rightarrow\rho^{0}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ decay. For comparison, we also list a subset of these couplings corresponding to the lowest $L$, used in the previous determinations [16, 17, 14] of the $X(3872)$ quantum numbers. | $B_{LS}$
---|---
$J^{PC}$ | Any $L$ value | Minimal $L$ value
$0^{-+}$ | $B_{11}$ | $B_{11}$
$0^{++}$ | $B_{00},B_{22}$ | $B_{00}$
$1^{-+}$ | $B_{10},B_{11},B_{12},B_{32}$ | $B_{10},B_{11},B_{12}$
$1^{++}$ | $B_{01},B_{21},B_{22}$ | $B_{01}$
$2^{-+}$ | $B_{11},B_{12},B_{31},B_{32}$ | $B_{11},B_{12}$
$2^{++}$ | $B_{02},B_{20},B_{21},B_{22},B_{42}$ | $B_{02}$
$3^{-+}$ | $B_{12},B_{30},B_{31},B_{32},B_{52}$ | $B_{12}$
$3^{++}$ | $B_{21},B_{22},B_{41},B_{42}$ | $B_{21},B_{22}$
$4^{-+}$ | $B_{31},B_{32},B_{51},B_{52}$ | $B_{31},B_{32}$
$4^{++}$ | $B_{22},B_{40},B_{41},B_{42},B_{62}$ | $B_{22}$

The function to be minimized is $-2\ln{\cal L}(J_{X},\alpha)\equiv-
s_{w}\,2\sum_{i=1}^{N_{\rm data}}w_{i}\ln{\cal P}(\Omega_{i}|J_{X},\alpha)$,
where ${\cal L}(J_{X},\alpha)$ is the unbinned likelihood, and $N_{\rm data}$
is the number of selected candidates. The background is subtracted using the
sPlot technique [32] by assigning a weight, $w_{i}$, to each candidate based
on its $\Delta M$ value (see Fig. 1). No correlations between $\Delta M$ and
$\Omega$ are observed. Prompt production of $X(3872)$ in $pp$ collisions gives
negligible contribution to the selected sample. Statistical fluctuations in
the background subtraction are taken into account in the log-likelihood value
via a constant scaling factor, $s_{w}=\sum_{i=1}^{N_{\rm
data}}w_{i}/\sum_{i=1}^{N_{\rm data}}{w_{i}}^{2}$. The efficiency
$\epsilon(\Omega)$ is not determined on an event-by-event basis, since it
cancels in the likelihood ratio except for the normalization integrals. A
large sample of simulated events, with uniform angular distributions, passed
through a full simulation of the detection and the data selection process, is
used to carry out the integration, $I(J_{X})\propto\sum_{i=1}^{N_{\rm
MC}}|{\cal M}(\Omega_{i}|J_{X})|^{2}$, where $N_{\rm MC}$ is the number of
reconstructed simulated events. The negative log-likelihood is minimized for
each $J_{X}$ value with respect to free $B_{LS}$ couplings, yielding their
estimated set of values $\hat{\alpha}$. Hereinafter, ${\cal
L}(J_{X})\equiv{\cal L}(J_{X},\hat{\alpha})$.

The $1^{++}$ hypothesis gives the highest likelihood value. From angular
momentum and parity conservation, there are two possible values of orbital
angular momentum in the $X(3872)$ decay for this $J^{PC}$ value, $L=0$ or $2$.
For the S–wave decay, the total spin of the $\rho^{0}$ and
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ mesons must be $S=1$; thus
$B_{01}$ is the only possible $LS$ amplitude. For the D–wave decay, two values
are possible, $S=1$ or $2$, corresponding to the amplitudes $B_{21}$ and
$B_{22}$, respectively. The squared magnitudes of both of these D–wave
amplitudes are consistent with zero, as demonstrated by the ratios
$|B_{21}|^{2}/|B_{01}|^{2}=0.002\pm 0.004$ and
$|B_{22}|^{2}/|B_{01}|^{2}=0.007\pm 0.008$. Overall, the D–wave significance
is only $0.8$ standard deviations as obtained by applying Wilks theorem to the
ratio of the likelihood values with the D–wave amplitudes floated in the fit
and with them fixed to zero. The total D–wave fraction depends on the $B_{LS}$
amplitudes, $f_{\rm D}\equiv\int\left|{\cal M}(\Omega)_{\rm
D}\right|^{2}d\Omega/\int\left|{\cal M}(\Omega)_{\rm S+D}\right|^{2}d\Omega$,
where ${\cal M}(\Omega)_{D}$ is the matrix element restricted to the $B_{21}$
and $B_{22}$ amplitudes only and ${\cal M}(\Omega)_{S+D}$ is the full matrix
element. To set an upper limit on $f_{\rm D}$, we populate the four-
dimensional space of complex $B_{21}$ and $B_{22}$ parameters with uniformly
distributed points in a large region around the $B_{21}$ and $B_{22}$ fit
values ($\pm 14$ standard deviations in each parameter). For each point we
determine the likelihood value from the data and an $f_{\rm D}$ value via
numerical integration of the matrix element squared. The distribution of
$f_{\rm D}$ values weighted by the likelihood values is shown in Fig. 2. It
peaks at $0.4\%$ with a non-Gaussian tail at higher values. An upper limit of
$f_{\rm D}<4\%$ at 95% C.L. is determined using a Bayesian approach.

Figure 2: Likelihood-weighted distribution of the D–wave fraction. The
distribution is normalized to unity.

The likelihood ratio $t\equiv-2\ln[{\cal L}(J_{X}^{\rm alt})/{\cal
L}(1^{++})]$ is used as a test variable to discriminate between the $1^{++}$
and alternative spin hypotheses considered ($J_{X}^{\rm alt}$). The values of
$t$ in the data ($t_{\rm data}$) are positive, favoring the $1^{++}$
assignment. They are incompatible with the distributions of $t$ observed in
experiments simulated under various $J_{X}^{\rm alt}$ hypotheses, as
illustrated in Fig. 3. To quantify these disagreements we calculate the
approximate significance of rejection (p-value) of $J_{X}^{\rm alt}$ as
$(t_{\rm data}-\langle t\rangle)/\sigma(t)$, where $\langle t\rangle$ and
$\sigma(t)$ are the mean and rms deviations of the $t$ distribution under the
$J_{X}^{\rm alt}$ hypothesis. In all spin configurations tested, we exclude
the alternative spin hypothesis with a significance of more than 16 standard
deviations. Values of $t$ in data are consistent with those expected in the
$1^{++}$ case as shown in Fig. 3, with fractions of simulated $1^{++}$
experiments with $t<t_{\rm data}$ in the 25%–91% range. Projections of the
data and of the fit ${\cal P}$ onto individual angles show good consistency
with the $1^{++}$ assignment as illustrated in Fig. 4. Inconsistency with the
other assignments is apparent when correlations between various angles are
exploited. For example, the data projection onto $\cos\theta_{X}$ is
consistent only with the $1^{++}$ fit projection after requiring
$|\cos\theta_{\rho}|>0.6$ (see Fig. 5), while inconsistency with the other
quantum number assignments is less clear without the $\cos\theta_{\rho}$
requirement.

Figure 3: Distributions of the test statistic $t\equiv-2\ln[{\cal
L}(J_{X}^{\rm alt}))/{\cal L}(1^{++})]$, for the simulated experiments under
the $J^{PC}=J_{X}^{\rm alt}$ hypothesis (blue solid histograms) and under the
$J^{PC}=1^{++}$ hypothesis (red dashed histograms). The values of the test
statistics for the data, $t_{\rm data}$, are shown by the solid vertical
lines.

The selection criteria are varied to probe for possible biases from the
background subtraction and the efficiency corrections. By requiring $Q<0.1$
GeV, the background level is reduced by more than a factor of two, while
losing only $20\%$ of the signal. By tightening the requirements on the
$p_{\rm T}$ of $\pi$, $K$ and $\mu$ candidates, we decrease the signal
efficiency by around 75% with similar reduction in the background level. In
all cases, the significance of the rejection of the disfavored hypotheses is
compatible with that expected from the simulation. Likewise, the best fit
$f_{\rm D}$ values determined for these subsamples of data change within the
expected statistical fluctuations and remain consistent with the upper limit
we have set.

Figure 4: Background-subtracted distributions of all angles for the data
(points with error bars) and for the $1^{++}$ fit projections (solid
histograms).  Figure 5: Background-subtracted distribution of $\cos\theta_{X}$
for candidates with $|\cos\theta_{\rho}|>0.6$ for the data (points with error
bars) compared to the expected distributions for various $X(3872)$ $J^{PC}$
assignments (solid histograms) with the $B_{LS}$ amplitudes obtained by the
fit to the data in the five-dimensional angular space. The fit displays are
normalized to the observed number of the signal events in the full angular
phase space.

In summary, the analysis of the angular correlations in $B^{+}\rightarrow
X(3872)K^{+}$,
$X(3872)\rightarrow\pi^{+}\pi^{-}{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}$, ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}\rightarrow\mu^{+}\mu^{-}$ decays, performed for the first time without
any assumption about the orbital angular momentum in the $X(3872)$ decay,
confirms that the eigenvalues of total angular momentum, parity and charge-
conjugation of the $X(3872)$ state are $1^{++}$. These quantum numbers are
consistent with those predicted by the molecular or tetraquark models and with
the $\chi_{c1}(2^{3}{\rm P}_{1})$ charmonium state [33], possibly mixed with a
molecule [10]. Other charmonium states are excluded. No significant D–wave
fraction is found, with an upper limit of $4\%$ at 95% C.L. The S–wave
dominance is expected in the charmonium or tetraquark models, in which the
$X(3872)$ state has a compact size. An extended size, as that predicted by the
molecular model, implies more favorable conditions for the D wave. However,
conclusive discrimination among models is difficult because quantitative
predictions are not available.

We express our gratitude to our colleagues in the CERN accelerator departments
for the excellent performance of the LHC. We thank the technical and
administrative staff at the LHCb institutes. We acknowledge support from CERN
and from the national agencies: CAPES, CNPq, FAPERJ and FINEP (Brazil); NSFC
(China); CNRS/IN2P3 (France); BMBF, DFG, HGF and MPG (Germany); INFN (Italy);
FOM and NWO (The Netherlands); MNiSW and NCN (Poland); MEN/IFA (Romania);
MinES and FANO (Russia); MinECo (Spain); SNSF and SER (Switzerland); NASU
(Ukraine); STFC (United Kingdom); NSF (USA). The Tier1 computing centres are
supported by IN2P3 (France), KIT and BMBF (Germany), INFN (Italy), NWO and
SURF (The Netherlands), PIC (Spain), GridPP (United Kingdom). We are indebted
to the communities behind the multiple open source software packages on which
we depend. We are also thankful for the computing resources and the access to
software R&D tools provided by Yandex LLC (Russia). Individual groups or
members have received support from EPLANET, Marie Skłodowska-Curie Actions and
ERC (European Union), Conseil général de Haute-Savoie, Labex ENIGMASS and
OCEVU, Région Auvergne (France), RFBR (Russia), XuntaGal and GENCAT (Spain),
Royal Society and Royal Commission for the Exhibition of 1851 (United
Kingdom).

## References

  * [1] Belle collaboration, S.-K. Choi et al., _Observation of a narrow charmonium-like state in exclusive $B^{\pm}\rightarrow K^{\pm}\pi^{+}\pi^{-}J/\psi$ decays_, Phys. Rev. Lett. 91 (2003) 262001, arXiv:hep-ex/0309032
  * [2] CDF collaboration, D. Acosta et al., _Observation of the narrow state $X(3872)\rightarrow J/\psi\pi^{+}\pi^{-}$ in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV_, Phys. Rev. Lett. 93 (2004) 072001, arXiv:hep-ex/0312021
  * [3] D0 collaboration, V. M. Abazov et al., _Observation and properties of the $X(3872)$ decaying to $J/\psi\pi^{+}\pi^{-}$ in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV_, Phys. Rev. Lett. 93 (2004) 162002, arXiv:hep-ex/0405004
  * [4] BaBar collaboration, B. Aubert et al., _Study of the $B^{-}\rightarrow J/\psi K^{-}\pi^{+}\pi^{-}$ decay and measurement of the $B^{-}\rightarrow X(3872)K^{-}$ branching fraction_, Phys. Rev. D71 (2005) 071103, arXiv:hep-ex/0406022
  * [5] LHCb collaboration, R. Aaij et al., _Observation of $X(3872)$ production in $pp$ collisions at $\sqrt{s}=7$ TeV_, Eur. Phys. J. C72 (2012) 1972, arXiv:1112.5310
  * [6] CMS collaboration, S. Chatrchyan et al., _Measurement of the $X(3872)$ production cross section via decays to ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}\pi\pi$ in $pp$ collisions at $\sqrt{(}s)=7$ TeV_, JHEP 04 (2013) 154, arXiv:1302.3968
  * [7] Particle Data Group, K. A. Olive et al., _Review of particle physics_ , Chin. Phys. C38 (2014) 090001, arXiv:1412.1408
  * [8] N. A. Tornqvist, _Isospin breaking of the narrow charmonium state of Belle at 3872-MeV as a deuson_ , Phys. Lett. B590 (2004) 209, arXiv:hep-ph/0402237
  * [9] L. Maiani, F. Piccinini, A. D. Polosa, and V. Riquer, _Diquark-antidiquarks with hidden or open charm and the nature of $X(3872)$_, Phys. Rev. D71 (2005) 014028, arXiv:hep-ph/0412098
  * [10] C. Hanhart, Y. S. Kalashnikova, and A. V. Nefediev, _Interplay of quark and meson degrees of freedom in a near-threshold resonance: multi-channel case_ , Eur. Phys. J. A47 (2011) 101, arXiv:1106.1185
  * [11] BaBar collaboration, B. Aubert et al., _Search for $B^{+}\rightarrow X(3872)K^{+}$, $X(3872)\rightarrow J/\psi\gamma$_, Phys. Rev. D74 (2006) 071101, arXiv:hep-ex/0607050
  * [12] Belle collaboration, V. Bhardwaj et al., _Observation of $X(3872)\rightarrow J/\psi\gamma$ and search for $X(3872)\rightarrow\psi^{\prime}\gamma$ in B decays_, Phys. Rev. Lett. 107 (2011) 091803, arXiv:1105.0177
  * [13] LHCb collaboration, R. Aaij et al., _Evidence for the decay $X(3872)\rightarrow\psi(2S)\gamma$_, Nucl. Phys. B886 (2014) 665, arXiv:1404.0275
  * [14] Belle collaboration, S.-K. Choi et al., _Bounds on the width, mass difference and other properties of $X(3872)\rightarrow\pi^{+}\pi^{-}J/\psi$ decays_, Phys. Rev. D84 (2011) 052004, arXiv:1107.0163
  * [15] CDF collaboration, A. Abulencia et al., _Measurement of the dipion mass spectrum in $X(3872)\rightarrow J/\psi\pi^{+}\pi^{-}$ decays._, Phys. Rev. Lett. 96 (2006) 102002, arXiv:hep-ex/0512074
  * [16] CDF collaboration, A. Abulencia et al., _Analysis of the quantum numbers $J^{PC}$ of the $X(3872)$_, Phys. Rev. Lett. 98 (2007) 132002, arXiv:hep-ex/0612053
  * [17] LHCb collaboration, R. Aaij et al., _Determination of the $X(3872)$ meson quantum numbers_, Phys. Rev. Lett. 110 (2013) 222001, arXiv:1302.6269
  * [18] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [19] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [20] T. Skwarnicki, A study of the radiative cascade transitions between the $\Upsilon^{\prime}$ and $\Upsilon$ resonances, PhD thesis, Institute of Nuclear Physics, Krakow, 1986, DESY-F31-86-02
  * [21] M. Jacob and G. C. Wick, _On the general theory of collisions for particles with spin_ , Annals Phys. 7 (1959) 404
  * [22] J. D. Richman, _An experimenter’s guide to the helicity formalism_ ,  1984, CALT-68-1148
  * [23] S. U. Chung, _General formulation of covariant helicity-coupling amplitudes_ , Phys. Rev. D57 (1998) 431
  * [24] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [25] I. Belyaev et al., _Handling of the generation of primary events in Gauss, the LHCb simulation framework_, Nuclear Science Symposium Conference Record (NSS/MIC) IEEE (2010) 1155
  * [26] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [27] GEANT4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [28] GEANT4 collaboration, S. Agostinelli et al., _GEANT4: A simulation toolkit_ , Nucl. Instrum. Meth. A506 (2003) 250
  * [29] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. of Phys. : Conf. Ser. 331 (2011) 032023
  * [30] LHCb collaboration, R. Aaij et al., _Observation of overlapping spin- $1$ and spin-$3$ $\bar{D}^{0}K^{-}$ resonances at mass $2.86$ GeV/$c^{2}$_, Phys. Rev. Lett. 113 (2014) 162001, arXiv:1407.7574
  * [31] LHCb collaboration, R. Aaij et al., _Dalitz plot analysis of $B^{0}_{s}\rightarrow\bar{D}^{0}K^{-}\pi^{+}$ decays_, Phys. Rev. D90 (2014) 072003, arXiv:1407.7712
  * [32] M. Pivk and F. R. Le Diberder, _sPlot: A statistical tool to unfold data distributions_ , Nucl. Instrum. Meth. A555 (2005) 356, arXiv:physics/0402083v3
  * [33] N. N. Achasov and E. V. Rogozina, _$X(3872)$ , $I^{G}(J^{PC})=0^{+}(1^{++})$, as the $\chi_{1c}(2P)$ charmonium_, arXiv:1501.03583

LHCb collaboration

R. Aaij38, B. Adeva37, M. Adinolfi46, A. Affolder52, Z. Ajaltouni5, S. Akar6,
J. Albrecht9, F. Alessio38, M. Alexander51, S. Ali41, G. Alkhazov30, P.
Alvarez Cartelle53, A.A. Alves Jr57, S. Amato2, S. Amerio22, Y. Amhis7, L.
An3, L. Anderlini17,g, J. Anderson40, M. Andreotti16,f, J.E. Andrews58, R.B.
Appleby54, O. Aquines Gutierrez10, F. Archilli38, P. d’Argent11, A.
Artamonov35, M. Artuso59, E. Aslanides6, G. Auriemma25,n, M. Baalouch5, S.
Bachmann11, J.J. Back48, A. Badalov36, C. Baesso60, W. Baldini16,38, R.J.
Barlow54, C. Barschel38, S. Barsuk7, W. Barter38, V. Batozskaya28, V.
Battista39, A. Bay39, L. Beaucourt4, J. Beddow51, F. Bedeschi23, I. Bediaga1,
L.J. Bel41, I. Belyaev31, E. Ben-Haim8, G. Bencivenni18, S. Benson38, J.
Benton46, A. Berezhnoy32, R. Bernet40, A. Bertolin22, M.-O. Bettler38, M. van
Beuzekom41, A. Bien11, S. Bifani45, T. Bird54, A. Birnkraut9, A. Bizzeti17,i,
T. Blake48, F. Blanc39, J. Blouw10, S. Blusk59, V. Bocci25, A. Bondar34, N.
Bondar30,38, W. Bonivento15, S. Borghi54, M. Borsato7, T.J.V. Bowcock52, E.
Bowen40, C. Bozzi16, S. Braun11, D. Brett54, M. Britsch10, T. Britton59, J.
Brodzicka54, N.H. Brook46, A. Bursche40, J. Buytaert38, S. Cadeddu15, R.
Calabrese16,f, M. Calvi20,k, M. Calvo Gomez36,p, P. Campana18, D. Campora
Perez38, L. Capriotti54, A. Carbone14,d, G. Carboni24,l, R. Cardinale19,j, A.
Cardini15, P. Carniti20, L. Carson50, K. Carvalho Akiba2,38, R. Casanova
Mohr36, G. Casse52, L. Cassina20,k, L. Castillo Garcia38, M. Cattaneo38, Ch.
Cauet9, G. Cavallero19, R. Cenci23,t, M. Charles8, Ph. Charpentier38, M.
Chefdeville4, S. Chen54, S.-F. Cheung55, N. Chiapolini40, M. Chrzaszcz40,26,
X. Cid Vidal38, G. Ciezarek41, P.E.L. Clarke50, M. Clemencic38, H.V. Cliff47,
J. Closier38, V. Coco38, J. Cogan6, E. Cogneras5, V. Cogoni15,e, L.
Cojocariu29, G. Collazuol22, P. Collins38, A. Comerma-Montells11, A.
Contu15,38, A. Cook46, M. Coombes46, S. Coquereau8, G. Corti38, M. Corvo16,f,
B. Couturier38, G.A. Cowan50, D.C. Craik48, A. Crocombe48, M. Cruz Torres60,
S. Cunliffe53, R. Currie53, C. D’Ambrosio38, J. Dalseno46, P.N.Y. David41, A.
Davis57, K. De Bruyn41, S. De Capua54, M. De Cian11, J.M. De Miranda1, L. De
Paula2, W. De Silva57, P. De Simone18, C.-T. Dean51, D. Decamp4, M.
Deckenhoff9, L. Del Buono8, N. Déléage4, D. Derkach55, O. Deschamps5, F.
Dettori38, B. Dey40, A. Di Canto38, F. Di Ruscio24, H. Dijkstra38, S.
Donleavy52, F. Dordei11, M. Dorigo39, A. Dosil Suárez37, D. Dossett48, A.
Dovbnya43, K. Dreimanis52, L. Dufour41, G. Dujany54, F. Dupertuis39, P.
Durante38, R. Dzhelyadin35, A. Dziurda26, A. Dzyuba30, S. Easo49,38, U.
Egede53, V. Egorychev31, S. Eidelman34, S. Eisenhardt50, U. Eitschberger9, R.
Ekelhof9, L. Eklund51, I. El Rifai5, Ch. Elsasser40, S. Ely59, S. Esen11, H.M.
Evans47, T. Evans55, A. Falabella14, C. Färber11, C. Farinelli41, N. Farley45,
S. Farry52, R. Fay52, D. Ferguson50, V. Fernandez Albor37, F. Ferrari14, F.
Ferreira Rodrigues1, M. Ferro-Luzzi38, S. Filippov33, M. Fiore16,38,f, M.
Fiorini16,f, M. Firlej27, C. Fitzpatrick39, T. Fiutowski27, P. Fol53, M.
Fontana10, F. Fontanelli19,j, R. Forty38, O. Francisco2, M. Frank38, C.
Frei38, M. Frosini17, J. Fu21, E. Furfaro24,l, A. Gallas Torreira37, D.
Galli14,d, S. Gallorini22,38, S. Gambetta19,j, M. Gandelman2, P. Gandini55, Y.
Gao3, J. García Pardiñas37, J. Garofoli59, J. Garra Tico47, L. Garrido36, D.
Gascon36, C. Gaspar38, U. Gastaldi16, R. Gauld55, L. Gavardi9, G. Gazzoni5, A.
Geraci21,v, D. Gerick11, E. Gersabeck11, M. Gersabeck54, T. Gershon48, Ph.
Ghez4, A. Gianelle22, S. Gianì39, V. Gibson47, L. Giubega29, V.V. Gligorov38,
C. Göbel60, D. Golubkov31, A. Golutvin53,31,38, A. Gomes1,a, C. Gotti20,k, M.
Grabalosa Gándara5, R. Graciani Diaz36, L.A. Granado Cardoso38, E. Graugés36,
E. Graverini40, G. Graziani17, A. Grecu29, E. Greening55, S. Gregson47, P.
Griffith45, L. Grillo11, O. Grünberg63, B. Gui59, E. Gushchin33, Yu. Guz35,38,
T. Gys38, C. Hadjivasiliou59, G. Haefeli39, C. Haen38, S.C. Haines47, S.
Hall53, B. Hamilton58, T. Hampson46, X. Han11, S. Hansmann-Menzemer11, N.
Harnew55, S.T. Harnew46, J. Harrison54, J. He38, T. Head39, V. Heijne41, K.
Hennessy52, P. Henrard5, L. Henry8, J.A. Hernando Morata37, E. van
Herwijnen38, M. Heß63, A. Hicheur2, D. Hill55, M. Hoballah5, C. Hombach54, W.
Hulsbergen41, T. Humair53, N. Hussain55, D. Hutchcroft52, D. Hynds51, M.
Idzik27, P. Ilten56, R. Jacobsson38, A. Jaeger11, J. Jalocha55, E. Jans41, A.
Jawahery58, F. Jing3, M. John55, D. Johnson38, C.R. Jones47, C. Joram38, B.
Jost38, N. Jurik59, S. Kandybei43, W. Kanso6, M. Karacson38, T.M. Karbach38,†,
S. Karodia51, M. Kelsey59, I.R. Kenyon45, M. Kenzie38, T. Ketel42, B.
Khanji20,38,k, C. Khurewathanakul39, S. Klaver54, K. Klimaszewski28, O.
Kochebina7, M. Kolpin11, I. Komarov39, R.F. Koopman42, P. Koppenburg41,38, M.
Korolev32, L. Kravchuk33, K. Kreplin11, M. Kreps48, G. Krocker11, P.
Krokovny34, F. Kruse9, W. Kucewicz26,o, M. Kucharczyk26, V. Kudryavtsev34, K.
Kurek28, T. Kvaratskheliya31, V.N. La Thi39, D. Lacarrere38, G. Lafferty54, A.
Lai15, D. Lambert50, R.W. Lambert42, G. Lanfranchi18, C. Langenbruch48, B.
Langhans38, T. Latham48, C. Lazzeroni45, R. Le Gac6, J. van Leerdam41, J.-P.
Lees4, R. Lefèvre5, A. Leflat32, J. Lefrançois7, O. Leroy6, T. Lesiak26, B.
Leverington11, Y. Li7, T. Likhomanenko65,64, M. Liles52, R. Lindner38, C.
Linn38, F. Lionetto40, B. Liu15, S. Lohn38, I. Longstaff51, J.H. Lopes2, P.
Lowdon40, D. Lucchesi22,r, H. Luo50, A. Lupato22, E. Luppi16,f, O. Lupton55,
F. Machefert7, F. Maciuc29, O. Maev30, K. Maguire54, S. Malde55, A. Malinin64,
G. Manca15,e, G. Mancinelli6, P. Manning59, A. Mapelli38, J. Maratas5, J.F.
Marchand4, U. Marconi14, C. Marin Benito36, P. Marino23,38,t, R. Märki39, J.
Marks11, G. Martellotti25, M. Martinelli39, D. Martinez Santos42, F. Martinez
Vidal66, D. Martins Tostes2, A. Massafferri1, R. Matev38, A. Mathad48, Z.
Mathe38, C. Matteuzzi20, A. Mauri40, B. Maurin39, A. Mazurov45, M. McCann53,
J. McCarthy45, A. McNab54, R. McNulty12, B. Meadows57, F. Meier9, M.
Meissner11, M. Merk41, D.A. Milanes62, M.-N. Minard4, D.S. Mitzel11, J. Molina
Rodriguez60, S. Monteil5, M. Morandin22, P. Morawski27, A. Mordà6, M.J.
Morello23,t, J. Moron27, A.B. Morris50, R. Mountain59, F. Muheim50, J.
Müller9, K. Müller40, V. Müller9, M. Mussini14, B. Muster39, P. Naik46, T.
Nakada39, R. Nandakumar49, I. Nasteva2, M. Needham50, N. Neri21, S. Neubert11,
N. Neufeld38, M. Neuner11, A.D. Nguyen39, T.D. Nguyen39, C. Nguyen-Mau39,q, V.
Niess5, R. Niet9, N. Nikitin32, T. Nikodem11, D. Ninci23, A. Novoselov35, D.P.
O’Hanlon48, A. Oblakowska-Mucha27, V. Obraztsov35, S. Ogilvy51, O.
Okhrimenko44, R. Oldeman15,e, C.J.G. Onderwater67, B. Osorio Rodrigues1, J.M.
Otalora Goicochea2, A. Otto38, P. Owen53, A. Oyanguren66, A. Palano13,c, F.
Palombo21,u, M. Palutan18, J. Panman38, A. Papanestis49, M. Pappagallo51, L.L.
Pappalardo16,f, C. Parkes54, G. Passaleva17, G.D. Patel52, M. Patel53, C.
Patrignani19,j, A. Pearce54,49, A. Pellegrino41, G. Penso25,m, M. Pepe
Altarelli38, S. Perazzini14,d, P. Perret5, L. Pescatore45, K. Petridis46, A.
Petrolini19,j, M. Petruzzo21, E. Picatoste Olloqui36, B. Pietrzyk4, T.
Pilař48, D. Pinci25, A. Pistone19, S. Playfer50, M. Plo Casasus37, T.
Poikela38, F. Polci8, A. Poluektov48,34, I. Polyakov31, E. Polycarpo2, A.
Popov35, D. Popov10, B. Popovici29, C. Potterat2, E. Price46, J.D. Price52, J.
Prisciandaro39, A. Pritchard52, C. Prouve46, V. Pugatch44, A. Puig Navarro39,
G. Punzi23,s, W. Qian4, R. Quagliani7,46, B. Rachwal26, J.H. Rademacker46, B.
Rakotomiaramanana39, M. Rama23, M.S. Rangel2, I. Raniuk43, N. Rauschmayr38, G.
Raven42, F. Redi53, S. Reichert54, M.M. Reid48, A.C. dos Reis1, S.
Ricciardi49, S. Richards46, M. Rihl38, K. Rinnert52, V. Rives Molina36, P.
Robbe7,38, A.B. Rodrigues1, E. Rodrigues54, J.A. Rodriguez Lopez62, P.
Rodriguez Perez54, S. Roiser38, V. Romanovsky35, A. Romero Vidal37, M.
Rotondo22, J. Rouvinet39, T. Ruf38, H. Ruiz36, P. Ruiz Valls66, J.J. Saborido
Silva37, N. Sagidova30, P. Sail51, B. Saitta15,e, V. Salustino Guimaraes2, C.
Sanchez Mayordomo66, B. Sanmartin Sedes37, R. Santacesaria25, C. Santamarina
Rios37, M. Santimaria18, E. Santovetti24,l, A. Sarti18,m, C. Satriano25,n, A.
Satta24, D.M. Saunders46, D. Savrina31,32, M. Schiller38, H. Schindler38, M.
Schlupp9, M. Schmelling10, T. Schmelzer9, B. Schmidt38, O. Schneider39, A.
Schopper38, M.-H. Schune7, R. Schwemmer38, B. Sciascia18, A. Sciubba25,m, A.
Semennikov31, I. Sepp53, N. Serra40, J. Serrano6, L. Sestini22, P. Seyfert11,
M. Shapkin35, I. Shapoval16,43,f, Y. Shcheglov30, T. Shears52, L. Shekhtman34,
V. Shevchenko64, A. Shires9, R. Silva Coutinho48, G. Simi22, M. Sirendi47, N.
Skidmore46, I. Skillicorn51, T. Skwarnicki59, E. Smith55,49, E. Smith53, J.
Smith47, M. Smith54, H. Snoek41, M.D. Sokoloff57,38, F.J.P. Soler51, F.
Soomro39, D. Souza46, B. Souza De Paula2, B. Spaan9, P. Spradlin51, S.
Sridharan38, F. Stagni38, M. Stahl11, S. Stahl38, O. Steinkamp40, O.
Stenyakin35, F. Sterpka59, S. Stevenson55, S. Stoica29, S. Stone59, B.
Storaci40, S. Stracka23,t, M. Straticiuc29, U. Straumann40, R. Stroili22, L.
Sun57, W. Sutcliffe53, K. Swientek27, S. Swientek9, V. Syropoulos42, M.
Szczekowski28, P. Szczypka39,38, T. Szumlak27, S. T’Jampens4, T. Tekampe9, M.
Teklishyn7, G. Tellarini16,f, F. Teubert38, C. Thomas55, E. Thomas38, J. van
Tilburg41, V. Tisserand4, M. Tobin39, J. Todd57, S. Tolk42, L. Tomassetti16,f,
D. Tonelli38, S. Topp-Joergensen55, N. Torr55, E. Tournefier4, S. Tourneur39,
K. Trabelsi39, M.T. Tran39, M. Tresch40, A. Trisovic38, A. Tsaregorodtsev6, P.
Tsopelas41, N. Tuning41,38, A. Ukleja28, A. Ustyuzhanin65,64, U. Uwer11, C.
Vacca15,e, V. Vagnoni14, G. Valenti14, A. Vallier7, R. Vazquez Gomez18, P.
Vazquez Regueiro37, C. Vázquez Sierra37, S. Vecchi16, J.J. Velthuis46, M.
Veltri17,h, G. Veneziano39, M. Vesterinen11, B. Viaud7, D. Vieira2, M. Vieites
Diaz37, X. Vilasis-Cardona36,p, A. Vollhardt40, D. Volyanskyy10, D. Voong46,
A. Vorobyev30, V. Vorobyev34, C. Voß63, J.A. de Vries41, R. Waldi63, C.
Wallace48, R. Wallace12, J. Walsh23, S. Wandernoth11, J. Wang59, D.R. Ward47,
N.K. Watson45, D. Websdale53, A. Weiden40, M. Whitehead48, D. Wiedner11, G.
Wilkinson55,38, M. Wilkinson59, M. Williams38, M.P. Williams45, M. Williams56,
F.F. Wilson49, J. Wimberley58, J. Wishahi9, W. Wislicki28, M. Witek26, G.
Wormser7, S.A. Wotton47, S. Wright47, K. Wyllie38, Y. Xie61, Z. Xu39, Z.
Yang3, X. Yuan34, O. Yushchenko35, M. Zangoli14, M. Zavertyaev10,b, L. Zhang3,
Y. Zhang3, A. Zhelezov11, A. Zhokhov31, L. Zhong3.

1Centro Brasileiro de Pesquisas Físicas (CBPF), Rio de Janeiro, Brazil
2Universidade Federal do Rio de Janeiro (UFRJ), Rio de Janeiro, Brazil
3Center for High Energy Physics, Tsinghua University, Beijing, China
4LAPP, Université Savoie Mont-Blanc, CNRS/IN2P3, Annecy-Le-Vieux, France
5Clermont Université, Université Blaise Pascal, CNRS/IN2P3, LPC, Clermont-
Ferrand, France
6CPPM, Aix-Marseille Université, CNRS/IN2P3, Marseille, France
7LAL, Université Paris-Sud, CNRS/IN2P3, Orsay, France
8LPNHE, Université Pierre et Marie Curie, Université Paris Diderot,
CNRS/IN2P3, Paris, France
9Fakultät Physik, Technische Universität Dortmund, Dortmund, Germany
10Max-Planck-Institut für Kernphysik (MPIK), Heidelberg, Germany
11Physikalisches Institut, Ruprecht-Karls-Universität Heidelberg, Heidelberg,
Germany
12School of Physics, University College Dublin, Dublin, Ireland
13Sezione INFN di Bari, Bari, Italy
14Sezione INFN di Bologna, Bologna, Italy
15Sezione INFN di Cagliari, Cagliari, Italy
16Sezione INFN di Ferrara, Ferrara, Italy
17Sezione INFN di Firenze, Firenze, Italy
18Laboratori Nazionali dell’INFN di Frascati, Frascati, Italy
19Sezione INFN di Genova, Genova, Italy
20Sezione INFN di Milano Bicocca, Milano, Italy
21Sezione INFN di Milano, Milano, Italy
22Sezione INFN di Padova, Padova, Italy
23Sezione INFN di Pisa, Pisa, Italy
24Sezione INFN di Roma Tor Vergata, Roma, Italy
25Sezione INFN di Roma La Sapienza, Roma, Italy
26Henryk Niewodniczanski Institute of Nuclear Physics Polish Academy of
Sciences, Kraków, Poland
27AGH - University of Science and Technology, Faculty of Physics and Applied
Computer Science, Kraków, Poland
28National Center for Nuclear Research (NCBJ), Warsaw, Poland
29Horia Hulubei National Institute of Physics and Nuclear Engineering,
Bucharest-Magurele, Romania
30Petersburg Nuclear Physics Institute (PNPI), Gatchina, Russia
31Institute of Theoretical and Experimental Physics (ITEP), Moscow, Russia
32Institute of Nuclear Physics, Moscow State University (SINP MSU), Moscow,
Russia
33Institute for Nuclear Research of the Russian Academy of Sciences (INR RAN),
Moscow, Russia
34Budker Institute of Nuclear Physics (SB RAS) and Novosibirsk State
University, Novosibirsk, Russia
35Institute for High Energy Physics (IHEP), Protvino, Russia
36Universitat de Barcelona, Barcelona, Spain
37Universidad de Santiago de Compostela, Santiago de Compostela, Spain
38European Organization for Nuclear Research (CERN), Geneva, Switzerland
39Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
40Physik-Institut, Universität Zürich, Zürich, Switzerland
41Nikhef National Institute for Subatomic Physics, Amsterdam, The Netherlands
42Nikhef National Institute for Subatomic Physics and VU University Amsterdam,
Amsterdam, The Netherlands
43NSC Kharkiv Institute of Physics and Technology (NSC KIPT), Kharkiv, Ukraine
44Institute for Nuclear Research of the National Academy of Sciences (KINR),
Kyiv, Ukraine
45University of Birmingham, Birmingham, United Kingdom
46H.H. Wills Physics Laboratory, University of Bristol, Bristol, United
Kingdom
47Cavendish Laboratory, University of Cambridge, Cambridge, United Kingdom
48Department of Physics, University of Warwick, Coventry, United Kingdom
49STFC Rutherford Appleton Laboratory, Didcot, United Kingdom
50School of Physics and Astronomy, University of Edinburgh, Edinburgh, United
Kingdom
51School of Physics and Astronomy, University of Glasgow, Glasgow, United
Kingdom
52Oliver Lodge Laboratory, University of Liverpool, Liverpool, United Kingdom
53Imperial College London, London, United Kingdom
54School of Physics and Astronomy, University of Manchester, Manchester,
United Kingdom
55Department of Physics, University of Oxford, Oxford, United Kingdom
56Massachusetts Institute of Technology, Cambridge, MA, United States
57University of Cincinnati, Cincinnati, OH, United States
58University of Maryland, College Park, MD, United States
59Syracuse University, Syracuse, NY, United States
60Pontifícia Universidade Católica do Rio de Janeiro (PUC-Rio), Rio de
Janeiro, Brazil, associated to 2
61Institute of Particle Physics, Central China Normal University, Wuhan,
Hubei, China, associated to 3
62Departamento de Fisica , Universidad Nacional de Colombia, Bogota, Colombia,
associated to 8
63Institut für Physik, Universität Rostock, Rostock, Germany, associated to 11
64National Research Centre Kurchatov Institute, Moscow, Russia, associated to
31
65Yandex School of Data Analysis, Moscow, Russia, associated to 31
66Instituto de Fisica Corpuscular (IFIC), Universitat de Valencia-CSIC,
Valencia, Spain, associated to 36
67Van Swinderen Institute, University of Groningen, Groningen, The
Netherlands, associated to 41

aUniversidade Federal do Triângulo Mineiro (UFTM), Uberaba-MG, Brazil
bP.N. Lebedev Physical Institute, Russian Academy of Science (LPI RAS),
Moscow, Russia
cUniversità di Bari, Bari, Italy
dUniversità di Bologna, Bologna, Italy
eUniversità di Cagliari, Cagliari, Italy
fUniversità di Ferrara, Ferrara, Italy
gUniversità di Firenze, Firenze, Italy
hUniversità di Urbino, Urbino, Italy
iUniversità di Modena e Reggio Emilia, Modena, Italy
jUniversità di Genova, Genova, Italy
kUniversità di Milano Bicocca, Milano, Italy
lUniversità di Roma Tor Vergata, Roma, Italy
mUniversità di Roma La Sapienza, Roma, Italy
nUniversità della Basilicata, Potenza, Italy
oAGH - University of Science and Technology, Faculty of Computer Science,
Electronics and Telecommunications, Kraków, Poland
pLIFAELS, La Salle, Universitat Ramon Llull, Barcelona, Spain
qHanoi University of Science, Hanoi, Viet Nam
rUniversità di Padova, Padova, Italy
sUniversità di Pisa, Pisa, Italy
tScuola Normale Superiore, Pisa, Italy
uUniversità degli Studi di Milano, Milano, Italy
vPolitecnico di Milano, Milano, Italy
†Deceased


# Scaling Advantage in Approximate Optimization with Quantum Annealing

Humberto Munoz Bauza  Department of Physics and Astronomy, University of
Southern California, Los Angeles, California 90089, USA  Center for Quantum
Information Science & Technology, University of Southern California, Los
Angeles, California 90089, USA  Daniel Lidar  Department of Physics and
Astronomy, University of Southern California, Los Angeles, California 90089,
USA  Center for Quantum Information Science & Technology, University of
Southern California, Los Angeles, California 90089, USA  Department of
Electrical Engineering, University of Southern California, Los Angeles,
California 90089, USA  Department of Chemistry, University of Southern
California, Los Angeles, California 90089, USA

###### Abstract

Quantum annealing is a heuristic optimization algorithm that exploits quantum
evolution to approximately find lowest energy states [1, 2]. Quantum annealers
have scaled up in recent years to tackle increasingly larger and more highly
connected discrete optimization and quantum simulation problems [3, 4, 5, 6,
7]. Nevertheless, despite numerous attempts, a computational quantum advantage
in exact optimization using quantum annealing hardware has so far remained
elusive [8, 9, 10, 11, 12, 13, 14, 15, 16]. Here, we present evidence for a
quantum annealing scaling advantage in approximate optimization. The advantage
is relative to the top classical heuristic algorithm: parallel tempering with
isoenergetic cluster moves (PT-ICM) [17]. The setting is a family of 2D spin-
glass problems with high-precision spin-spin interactions. To achieve this
advantage, we implement quantum annealing correction (QAC) [18]: an embedding
of a bit-flip error-correcting code with energy penalties that leverages the
properties of the D-Wave Advantage quantum annealer to yield over $1,300$
error-suppressed logical qubits on a degree-$5$ interaction graph. We generate
random spin-glass instances on this graph and benchmark their time-to-epsilon,
a generalization of the time-to-solution metric [8] for low-energy states. We
demonstrate that with QAC, quantum annealing exhibits a scaling advantage over
PT-ICM at sampling low energy states with an optimality gap [19] of at least
$1.0\%$. This amounts to the first demonstration of an algorithmic quantum
speedup in approximate optimization.

The pursuit of a quantum speedup using quantum processors is a major theme in
modern physics and computer science. Two of the leading application areas are
discrete optimization and quantum simulation. In the latter context,
impressive recent quantum annealing (QA) advances have been made for fast,
coherent anneals lasting on the order of the single superconducting flux qubit
coherence time [20, 21]. While this diabatic approach is considered promIsing
[22], it cannot be expected to scale up without the introduction of error
correction or suppression, as decoherence and control errors pose
insurmountable challenges to Hamiltonian quantum computation models [23, 24,
25, 26, 27, 28], just as they do for gate-model quantum computers. In the
absence of a fault-tolerance threshold theorem [29] for QA, a variety of
Hamiltonian error suppression techniques have been proposed and analyzed as
ways to reduce the error rates of this computational model and the closely
related model of adiabatic quantum computation [30, 31, 32, 33, 34, 35, 36],
providing tools towards scalability.

However, despite these theoretical advances, there are currently practical
limitations to the types and locality of programmable interactions in the
Hamiltonians of quantum annealing hardware, even as new devices have started
to emerge [37, 38, 39, 16]. To address these limitations, quantum annealing
correction (QAC) [18] was developed as a realizable error suppression method
for quantum annealing, targeting the available and restricted set of control
operations in quantum annealers. QAC has been demonstrated to enhance the
success probability of quantum annealing and mitigate the analog control
problem [18, 40, 41]. The QAC method is based on a repetition-code encoding of
the Hamiltonian and does not fully realize a Hamiltonian stabilizer code.
Despite this, it has been shown theoretically to increase the energy gap of
the encoded Hamiltonian and reduce tunneling barriers, thus softening the
onset of the associated critical dynamics as well as lowering the effective
temperature [42].

Here, departing from the traditional paradigm of using QA for exact
optimization, we demonstrate—by incorporating QAC—the first genuine scaling
advantage in _approximate_ optimization (low-energy sampling) using a quantum
annealer. Even approximate optimization can be computationally hard unless
$\mathrm{P=NP}$ [43, 44], so we do not expect the advantage we exhibit to
amount to more than a polynomial speedup. However, whereas the scaling
advantages reported in previous work were relative to simulated annealing [13,
16], the advantage we find here is over the best currently available general
heuristic classical optimization method: parallel tempering with isoenergetic
cluster moves (PT-ICM) [17]. This result is enabled by implementing QAC on the
D-Wave Advantage quantum annealer for the Sidon-set spin glass instance class
[9], embedded on the logical graph formed after the encoding step. The
advantage of quantum annealing over PT-ICM is diminished without QAC, thus
highlighting the crucial role of quantum error suppression.

Figure 1:  Embedding of QAC on the Pegasus hardware graph. (a) Close-up of the
embedding of three clusters of two coupled logical qubits each. Physical
qubits are denoted by circles and couplings by lines. Penalty qubits are
marked by a P and the three penalty couplings (thin lines) to their
corresponding physical data qubits. Thick lines indicate the physical
couplings between the corresponding physical qubits of different logical
qubits. Only a subset of all possible couplings are shown. (b) A zoomed-out
view from (a) showing the logical qubits (circles) and the logical couplings
induced by the QAC embedding. Thick lines indicate the logical couplings shown
in (a), whose type is colored by direction: horizontal/vertical/diagonal (long
or short). (c) A zoomed-out view from (b) showing a $4\times 4$ induced
logical graph of QAC. The logical graph is equivalent to a honeycomb graph
with additional non-planar bonds. The induced logical graph of the D-Wave
Advantage 4.1 QPU has a maximum of $1322$ logical qubits; the largest
available side length is $L=15$.

Quantum annealing
The D-Wave quantum processing unit (QPU) uses superconducting flux qubits to
implement the transverse field Ising Hamiltonian

$H(s)=-A(s)\sum_{i\in\mathcal{V}}\sigma^{x}_{i}+B(s)H_{z},$ (1)

where $\mathcal{V}$ is the vertex set of the hardware graph of the QPU, $i$ is
the qubit index, $\sigma^{x}_{i}$ are Pauli matrices, and $A(s)$ and $B(s)$
are the annealing schedules, respectively decreasing to $0$ and increasing
from $0$ with $s:0\to 1$. $H_{z}$ is the Ising problem Hamiltonian:

$H_{z}=\sum_{i\in\mathcal{V}}h_{i}\sigma_{i}^{z}+\sum_{\\{i,j\\}\in\mathcal{E}}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z},$
(2)

where $h_{i}$ and $J_{ij}$ are programmable local fields and couplings,
respectively, and $\mathcal{E}$ is the edge set of the hardware graph. Many
NP-complete and NP-hard problems can be mapped to $H_{z}$ [45] by minor-
embedding onto the hardware graph. We performed QA experiments on the D-Wave
Advantage 4.1 QPU accessed through the D-Wave Leap cloud interface, featuring
the Pegasus hardware graph [46].

Quantum annealing correction
We implement the $[[3,1,3]]_{1}$ QAC encoding introduced in Ref. [18], which
encodes a logical qubit into three physical “data qubits,” each of which is
coupled to the same additional “energy penalty qubit” with a fixed coupling
strength $J_{p}$; the logical qubit is decoded via a majority vote on the data
qubits. The logical subgraph induced by the QAC encoding on the Pegasus graph
has a bulk degree of $5$ and admits native loops of length $5$. Fig. 1
illustrates the encoding and the induced logical graph. All previous QAC
results were obtained using the “Chimera” hardware graph of the previous
generation of D-Wave QPUs, which has degree $3$ and no odd-length loops. The
features of the induced Pegasus logical graph permit the benchmarking of
significantly harder spin-glass instances than was possible on Chimera. The
induced logical graphs we examine have side length $L\in[5,15]$, corresponding
to a problem size range of $N\in[142,1322]$ logical, error-corrected qubits.

Problems on the logical QAC graph can also be encoded directly by setting
$J_{p}=0$, resulting in three uncoupled and unprotected, parallel classical
copies of the problem instance. We then extract the energies of all three
copies as independent annealing samples and denote this “unprotected” QA
method by U3. We use the U3 method below to test whether QAC has a genuine
advantage over simple classical repetition coding.

Spin-glass instances
We generate random native spin-glass instances on the induced logical graph.
These types of instances have been widely used to benchmark the previous
D-Wave QPUs (with the Chimera hardware graph) against classical algorithms [4,
8, 9]. We tested three types of spin-glass disorder: binomial, where $J_{ij}$
was randomly selected as $\pm 1$ with equal probability, Sidon-28 (S28) [9],
where $J_{ij}$ was randomly sampled from the set
$\pm\\{8/28,13/28,19/28,1\\}$, and finally range 6 (R6) disorder, where
$J_{ij}\in\pm\\{1/6,\ldots,6/6\\}$. In a Sidon set, the sum of any two set
members gives a number that is not part of the set. Moreover, no five numbers
from the S28 set sum to zero, which prevents the occurrence of “floppy” qubits
[47] given the bulk degree-$5$ of the Pegasus graph. Binomial disorder
generally admits a degenerate ground state, simplifying the optimization
problem. In contrast, the S28 disorder can yield instances with a unique
ground state [26]. The ground states are robust to small errors in the
implementation of the $J_{ij}$ values in the binomial disorder case, but this
is not the case when high precision in implementing the $J_{ij}$ values is
required (as for Sidon disorder). The latter case is expected to benefit more
from QAC than the former [41]. From here on, we focus on the S28 case; see
Methods.

Figure 2:  Time-to-epsilon scaling for QAC, U3, and PT-ICM for Sidon-28 (S28)
spin-glass disorder. The bottom panels show the optimal annealing times of U3
and QAC. The top panels show the $\operatorname{\text{TT}\varepsilon}$ results
for the corresponding optimal annealing times, along with optimized PT-ICM
results. The straight lines are best fits assuming power law scaling
$\operatorname{\text{TT}\varepsilon}=cN^{\alpha}$ and the accompanying numbers
are the corresponding slopes $\alpha$. As indicated in the legends, the target
$\varepsilon$ increases from left to right, with a corresponding improvement
in quantum annealing’s performance: PT-ICM is better at low $\varepsilon$, but
for $\varepsilon=1\%$, while U3 is still worse than PT-ICM, QAC already
outperforms PT-ICM. At $\varepsilon=1.25\%$, both QAC and U3 have better
scaling than PT-ICM. For higher residual energy targets, the scaling of QAC
becomes unreliable since we can no longer guarantee that the optimal $t_{f}$
has been identified. We used $t_{f}\in[0.5,27]\,\mathrm{\mu s}$ and
$N\in[142,1322]$. Error bars for $\operatorname{\text{TT}\varepsilon}$ data
points are twice the standard error of the parameter estimate calculated using
bootstrapping. Filled (open) circles correspond to a $P$ value of $0.05$
($0.20$) that $t_{f}^{\text{opt}}(N)>t_{f}^{\min}$ for the corresponding
problem size $N$ (see Methods for details). We use only those $t_{f}$ values
for which $P\leq 0.05$ to compute the slopes.

Time-to-epsilon metric
The standard performance metric for exact optimization using heuristic solvers
is the time-to-solution (TTS), which quantifies the time to reach the ground
state at least once with $99\%$ probability, given that the success
probability per run is $p$: $\text{TTS}=t_{f}R$, where each annealing run
lasts time $t_{f}$ and $R=\frac{\log(1-0.99)}{\log(1-p)}$ is the expected
number of runs [8]. To address approximate optimization, we instead consider
the time required to reach an energy within a fraction $\varepsilon$ of the
ground state energy $E_{0}$, and define the time-to-epsilon for an instance as

$\operatorname{\text{TT}\varepsilon}=t_{f}\frac{\log(1-0.99)}{\log(1-p_{E\leq
E_{0}+\varepsilon\absolutevalue{E_{0}}})},$ (3)

where $p_{E\leq E_{0}+\varepsilon\absolutevalue{E_{0}}}$ is the probability
that the energy $E$ of a sample is no more than
$\varepsilon\absolutevalue{E_{0}}$ above $E_{0}$. The TTS is the special case
$\varepsilon=0$. In a mixed integer programming optimization context,
$\varepsilon$ is known as the _optimality gap_ [19], which is how we refer to
it here. The ground state energies are known for our instances (see Methods ),
so $\varepsilon$ is exactly calculated for each sample rather than bounded. An
alternative time metric is the residual energy density from the ground state
[21]; we focus on the optimality gap due to its relevance to benchmarking
approximate optimization algorithms.

We define $[\operatorname{\text{TT}\varepsilon}]_{q}$ of an instance class as
the $q$-th quantile of $\operatorname{\text{TT}\varepsilon}$ over the entire
instance class. Here, we focus only on the median quantile, $q=0.5$, denoted
$[\operatorname{\text{TT}\varepsilon}]_{\mathrm{Med}}$. For a given disorder,
instance size, and $\varepsilon$-target, we find the annealing time $t_{f}$
(and penalty strength for QAC) that minimizes
$[\operatorname{\text{TT}\varepsilon}]_{\mathrm{Med}}$. We restrict the
penalty coupling strengths to the set $J_{p}\in\\{0.1,0.2,0.3\\}$ to reduce
resource requirements for parameter optimization, as $J_{p}=0.2$ is the
penalty strength that most frequently optimizes the success probabilities of
individual instances, and the dependence on $J_{p}$ above $0.2$ becomes weak.

Fitting the $\operatorname{\text{TT}\varepsilon}$
Below, we fit the $\operatorname{\text{TT}\varepsilon}$ to a power law:
$\operatorname{\text{TT}\varepsilon}(N)=cN^{\alpha}$, where $\alpha$ is the
_scaling exponent_ , the quantity we use to quantify the scaling of the
different algorithms we compare. The choice of a power law is motivated by the
existence of an $O(N)$ classical algorithm for the residual energy density; we
describe such an algorithm in Methods (though this algorithm is utterly
impractical due to its huge prefactor). Due to the power law fit, we should
account for factors that can modify the scaling exponent. Indeed, we could use
all $N_{\max}$ qubits of the QPU and embed ${N_{\max}}/{N}$ parallel copies of
each problem of size $N$, then select the best of these copies. Since, in
reality, we work with only one copy due to a small fraction of the qubits and
couplers being absent, we multiply the $\operatorname{\text{TT}\varepsilon}$
by a factor of ${N}/{N_{\max}}$ [8]. The U3
$\operatorname{\text{TT}\varepsilon}$ is similarly multiplied by a constant
factor of $3/4$ since, due to a lack of needed couplings, each instance is
repeated only over the data qubits, thus leaving $1/4$ of the available
(penalty) qubits unused.

Parallel tempering algorithm
Our baseline classical algorithm is parallel tempering with isoenergetic
cluster moves (PT-ICM) [17]. The runtime of this algorithm has the best
scaling with problem size known in the task of finding the ground state of
various benchmark problems on D-Wave QPUs [12, 14], with the only known
exception being certain XORSAT instances for which highly specialized solvers
have been developed [48, 15]. Our optimization of the algorithmic parameters
of PT-ICM is described in Methods.

Results
It is well known that the TTS metric generates unreliable results unless the
annealing time is optimized for each size $N$ [8, 13]. This is because an
artificially high TTS at small $N$ results in an overly flat TTS scaling. The
same considerations apply to the $\operatorname{\text{TT}\varepsilon}$, so
here we find the annealing time $t_{f}$ that minimizes
$[\operatorname{\text{TT}\varepsilon}]_{\mathrm{Med}}$ for each $N$—denoted
$t_{f}^{\text{opt}}(N)$—and report the resulting median
$\operatorname{\text{TT}\varepsilon}$ and its scaling estimate for QAC and U3
in Fig. 2, along with the analogously optimized PT-ICM results.

The shortest available annealing time on the D-Wave Advantage QPU accessed via
Leap is $t_{f}^{\min}=0.5\,\mathrm{\mu s}$, and the bottom panels of Fig. 2
show that as the target residual energy density is increased, progressively
larger problem sizes are needed to ensure that
$t_{f}^{\text{opt}}(N)>t_{f}^{\min}$. We cannot rule out that with access to
lower annealing times, one would find $t_{f}^{\text{opt}}(N)<t_{f}^{\min}$ for
all $N$ values where we empirically find $t_{f}^{\text{opt}}(N)=t_{f}^{\min}$.
We thus formulate a null hypothesis for each $N$ that
$t_{f}^{\text{opt}}(N)\leq t_{f}^{\min}$ and compute a $P$ value as the
empirical number of bootstrap samples whose
$t_{f}^{\text{opt}}(N)=t_{f}^{\min}$, out of a total of $200$ samples (see
Methods for details of our statistical analysis). To compute the
$\operatorname{\text{TT}\varepsilon}$ scaling, i.e., the slope $\alpha$ in a
fit to $\operatorname{\text{TT}\varepsilon}=cN^{\alpha}$, for each
$\varepsilon$ we use only those $t_{f}^{\text{opt}}(N)$ values whose $P<0.05$
(filled circles in Fig. 2). We can thus be confident that the reported slopes
reflect the true scaling of U3 and QAC.

Our first observation is that _the QAC scaling is always better than the U3
scaling_ , which is consistent with previous studies concerning the effect of
analog coupling errors (“$J$-chaos”) on the TTS for spin-glass instances [41].
Such errors are expected for the S28 instances due to the relatively high
precision their specification requires.

Second, we observe that U3 and QAC reduce the absolute algorithmic runtime by
four orders of magnitude compared to PT-ICM. However, this is not a scaling
advantage, and since our PT-ICM calculations could be sped up by employing
faster classical processors, we do not consider this a robust finding.
Similarly, we exclude the programming and readout time used by the D-Wave QPU.
The primary source of overhead is the readout time per sample, which scales
with problem size and can reach $200\,\mathrm{\mu s}$ per sample for this QPU.
This timing varies by hardware generation, and its inclusion obfuscates the
dominant scaling source.

Third, and most significantly, we observe that _as the target optimality gap
increases QA’s scaling overtakes PT-ICM_. Notably, at $\varepsilon=1\%$, QAC
exhibits a scaling exponent of $1.69\pm 0.12$, compared to PT-ICM’s $1.93\pm
0.03$, and at $\varepsilon=1.25\%$, QAC and U3 exhibit scaling exponents of
$1.15\pm 0.22$ and $1.76\pm 0.06$, respectively, compared to PT-ICM’s $1.87\pm
0.02$. The scaling of U3 at optimal annealing times continues to decrease as
$\varepsilon$ increases; at $\varepsilon=1.5\%$ (not shown), the scaling
exponents for U3 and PT-ICM become $1.60\pm 0.07$ and $1.86\pm 0.04$,
respectively. _This is robust evidence of a quantum annealing scaling
advantage over the best available classical heuristic optimization algorithm_.

We are unable to determine the scaling of QAC for $\varepsilon>1.25\%$, as we
cannot confirm that $t_{f}^{\text{opt}}(N)>t_{f}^{\min}$ for any $N$; as can
be seen in Fig. 2, already for $\varepsilon=1.25\%$ only the largest three $N$
values satisfy the $P<0.05$ criterion. However, given the consistently better
scaling of QAC for lower values of $\varepsilon$, where
$t_{f}^{\text{opt}}(N)>t_{f}^{\min}$ for QAC over a range of problem sizes, it
is reasonable to conclude that the QAC scaling would be a further improvement
over U3 if its true $t_{f}^{\text{opt}}(N)$ could be established for
$\varepsilon>1.25\%$; this would require access to shorter annealing times or
larger system sizes.

We note that it is unsurprIsing that, given a sufficiently large target
optimality gap, the D-Wave QPU returns sample energies within that gap.
Similarly, we can expect that for large enough $\varepsilon$, even simulated
annealing or greedy descent will be nearly guaranteed success in polynomial
time. The significance of our result is that _QA reaches near-linear scaling
at a smaller optimality gap target than PT-ICM_. Thus, we refer to this result
as an _approximate optimization advantage for quantum annealing_. We also note
that Ref. [21] similarly reported a QA optimization advantage for the residual
energy density (for 3D spin glasses), but this was done at a _fixed_ problem
size (of $N=5374$ physical qubits), and was instead concerned with the
convergence of the residual energy density with the annealing time. We
reemphasize that, in contrast, we are reporting a scaling advantage as a
function of problem size, the proper context for quantum speedup claims. We
explain in Methods that our results stand regardless of whether an optimality
gap or residual energy density target is chosen for benchmarking.

Dynamical critical scaling
As an additional perspective on the different quantum annealing dynamics
resulting from error suppression, we study the difference in the dynamical
critical scaling under the Kibble-Zurek (KZ) scaling ansatz [49, 50], where an
annealing time of $t_{f}\sim L^{\mu}$ is required to suppress diabatic
excitations with a correlation length of $L$, and where $\mu$ is the dynamical
critical exponent or KZ exponent. We examined signatures of dynamical critical
scaling by calculating the Binder cumulant
$U=\frac{1}{2}\left(3-{\expectationvalue{q^{4}}}/{\expectationvalue{q^{2}}^{2}}\right)$,
where $\expectationvalue{\cdot}$ denotes the sample average, either from U3 or
after decoding with QAC, and $U$ is averaged over all instances for each
system size $N$ and annealing time $t_{f}$. Here
$q=\frac{1}{N}\sum_{i=1}^{N}\sigma^{z}_{i}\sigma^{z\prime}_{i}$ is the overlap
between two replicas, i.e., independently annealed $N$-spin states
$\\{\sigma^{z}_{i}\\}$ and $\\{\sigma^{z\prime}_{i}\\}$ of a given disorder
realization (set of couplings $J_{ij}$).

We performed a finite-size scaling analysis and data collapse, and the
collapsed Binder cumulant for the S28 instances is shown in Fig. 3. We find
that $\mu_{\text{QAC}}=4.81\pm 0.22$ (at a penalty strength of $0.1$) compared
to $\mu_{\text{U3}}=7.53\pm 0.47$. The reduction of $\mu$ is lost when
$\lambda=0.2$ (not shown), suggesting that $\lambda=0.1$ is optimal in the
sense of diabatic error suppression.

This effect indicates that QAC is much more effective at suppressing diabatic
excitations. I.e., at equal annealing times, the dynamics are more adiabatic
under QAC, in agreement with theory [42]. The significant reduction in $\mu$
suggests that in addition to $J$-chaos suppression, diabatic error suppression
by QAC is responsible for the improved $\operatorname{\text{TT}\varepsilon}$
and shorter optimal annealing times. Additional context is provided in
Methods.

Figure 3:  Collapsed Binder cumulant for samples collected for U3 and QAC with
a penalty strength of $\lambda=0.1$. Here $t_{f}\in[0.5,5.2]\mu$s and
$L\in[5,15]$.

Discussion
Using the largest available quantum annealer to date we have demonstrated an
approximate optimization time-scaling advantage for quantum annealing on a
family of spin-glass problems with low ground-state degeneracy and high-
precision couplings. Our demonstration involves up to $1322$ logical qubits,
the largest number to date in an error-corrected setting. Our key result is
the demonstration of an _algorithmic_ quantum scaling advantage. The advantage
is relative to PT-ICM, the best classical heuristic algorithm currently known
for such spin glass problems, and appears at optimality gaps $\gtrsim 1\%$.

There are a few limitations to the scope of our conclusion. First, our results
do not imply that finding states within small, constant gaps (and indeed
finding the ground state itself) is easy for quantum annealing, nor do they
imply that all spin glass problems are amenable to an approximate optimization
scaling advantage via quantum annealing. Second, being finite-range and two-
dimensional limits the range of applications of the problem family we have
studied here. To achieve an algorithmic quantum advantage in an application
setting, the next challenge for quantum optimization is demonstrating a
hardware-scalable advantage in densely connected problems at sufficiently
small optimality gaps.

Acknowledgements
We thank Dr. Evgeny Mozgunov for suggesting the use of the time-to-epsilon
metric and Dr. Victor Kasatkin for discussions about the theoretical lower
bound of finite-range approximate optimization. We also thank Dr. Mohammad
Amin, Dr. Carleton Coffrin, Dr. Itay Hen, and Dr. Tameem Albash for various
helpful discussions and suggestions. This material is based upon work
supported by the Defense Advanced Research Projects Agency (DARPA) under
Agreement No. HR00112190071. This research was also supported by the ARO MURI
grant W911NF-22-S-0007. The authors acknowledge the Center for Advanced
Research Computing (CARC) at the University of Southern California for
computing resources.

## References

  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse Ising model, Physical Review E 58, 5355 (1998).
  * Hauke  _et al._ [2020] P. Hauke, H. G. Katzgraber, W. Lechner, H. Nishimori, and W. D. Oliver, Perspectives of quantum annealing: Methods and implementations, Reports on Progress in Physics (2020).
  * Johnson  _et al._ [2011] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, Quantum annealing with manufactured spins, Nature 473, 194 (2011).
  * Boixo  _et al._ [2014] S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Evidence for quantum annealing with more than one hundred qubits, Nature Physics 10, 218 (2014).
  * King  _et al._ [2018] A. D. King, J. Carrasquilla, J. Raymond, I. Ozfidan, E. Andriyash, A. Berkley, M. Reis, T. Lanting, R. Harris, F. Altomare, K. Boothby, P. I. Bunyk, C. Enderud, A. Fréchette, E. Hoskinson, N. Ladizinsky, T. Oh, G. Poulin-Lamarre, C. Rich, Y. Sato, A. Y. Smirnov, L. J. Swenson, M. H. Volkmann, J. Whittaker, J. Yao, E. Ladizinsky, M. W. Johnson, J. Hilton, and M. H. Amin, Observation of topological phenomena in a programmable lattice of 1,800 qubits, Nature 560, 456 (2018).
  * Harris  _et al._ [2018] R. Harris, Y. Sato, A. J. Berkley, M. Reis, F. Altomare, M. H. Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, S. Huang, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, T. Lanting, R. Li, T. Medina, R. Molavi, R. Neufeld, T. Oh, I. Pavlov, I. Perminov, G. Poulin-Lamarre, C. Rich, A. Smirnov, L. Swenson, N. Tsai, M. Volkmann, J. Whittaker, and J. Yao, Phase transitions in a programmable quantum spin glass simulator, Science 361, 162 (2018).
  * Weinberg  _et al._ [2020] P. Weinberg, M. Tylutki, J. M. Rönkkö, J. Westerholm, J. A. Åström, P. Manninen, P. Törmä, and A. W. Sandvik, Scaling and diabatic effects in quantum annealing with a d-wave device, Physical Review Letters 124, 090502 (2020).
  * Ronnow  _et al._ [2014] T. F. Ronnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Martinis, D. A. Lidar, and M. Troyer, Defining and detecting quantum speedup, Science 345, 420 (2014).
  * Katzgraber  _et al._ [2015] H. G. Katzgraber, F. Hamze, Z. Zhu, A. J. Ochoa, and H. Munoz-Bauza, Seeking Quantum Speedup Through Spin Glasses: The Good, the Bad, and the Ugly, Physical Review X 5, 031026 (2015).
  * Venturelli  _et al._ [2015] D. Venturelli, S. Mandrà, S. Knysh, B. O’Gorman, R. Biswas, and V. Smelyanskiy, Quantum optimization of fully connected spin glasses, Phys. Rev. X 5, 031040 (2015).
  * Denchev  _et al._ [2016] V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Babbush, V. Smelyanskiy, J. Martinis, and H. Neven, What is the computational value of finite-range tunneling?, Phys. Rev. X 6, 031015 (2016).
  * Mandrà  _et al._ [2017] S. Mandrà, H. G. Katzgraber, and C. Thomas, The pitfalls of planar spin-glass benchmarks: RaIsing the bar for quantum annealers (again), Quantum Science and Technology 2, 038501 (2017).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Demonstration of a Scaling Advantage for a Quantum Annealer over Simulated Annealing, Physical Review X 8, 031016 (2018).
  * Mandrà and Katzgraber [2018] S. Mandrà and H. G. Katzgraber, A deceptive step towards quantum speedup detection, Quantum Science and Technology 3, 04LT01 (2018).
  * Kowalsky  _et al._ [2022] M. Kowalsky, T. Albash, I. Hen, and D. A. Lidar, 3-regular three-xorsat planted solutions benchmark of classical and quantum heuristic optimizers, Quantum Science and Technology 7, 025008 (2022).
  * Ebadi  _et al._ [2022] S. Ebadi, A. Keesling, M. Cain, T. T. Wang, H. Levine, D. Bluvstein, G. Semeghini, A. Omran, J. Liu, R. Samajdar, X.-Z. Luo, B. Nash, X. Gao, B. Barak, E. Farhi, S. Sachdev, N. Gemelke, L. Zhou, S. Choi, H. Pichler, S. Wang, M. Greiner, V. Vuletic, and M. D. Lukin, Quantum Optimization of Maximum Independent Set using Rydberg Atom Arrays, Science 376, 1209 (2022).
  * Zhu  _et al._ [2015] Z. Zhu, A. J. Ochoa, and H. G. Katzgraber, Efficient Cluster Algorithm for Spin Glasses in Any Space Dimension, Physical Review Letters 115, 077201 (2015).
  * Pudenz  _et al._ [2014] K. L. Pudenz, T. Albash, and D. A. Lidar, Error-corrected quantum annealing with hundreds of qubits, Nature Communications 5, 3243 (2014).
  * [19] Gurobi Optimization: MIPGap.
  * King  _et al._ [2022] A. D. King, S. Suzuki, J. Raymond, A. Zucca, T. Lanting, F. Altomare, A. J. Berkley, S. Ejtemaee, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, J. D. Whittaker, J. Yao, R. Harris, D. A. Lidar, H. Nishimori, and M. H. Amin, Coherent quantum annealing in a programmable 2,000 qubit Ising chain, Nature Physics 18, 1324 (2022).
  * King  _et al._ [2023] A. D. King, J. Raymond, T. Lanting, R. Harris, A. Zucca, F. Altomare, A. J. Berkley, K. Boothby, S. Ejtemaee, C. Enderud, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, R. Molavi, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, N. Tsai, M. Volkmann, J. D. Whittaker, J. Yao, A. W. Sandvik, and M. H. Amin, Quantum critical dynamics in a 5,000-qubit programmable spin glass, Nature 617, 61 (2023).
  * Crosson and Lidar [2021] E. J. Crosson and D. A. Lidar, Prospects for quantum enhancement with diabatic quantum annealing, Nature Reviews Physics 3, 466 (2021).
  * Childs  _et al._ [2001] A. M. Childs, E. Farhi, and J. Preskill, Robustness of adiabatic quantum computation, Phys. Rev. A 65, 012322 (2001).
  * Amin  _et al._ [2009] M. H. S. Amin, D. V. Averin, and J. A. Nesteroff, Decoherence in adiabatic quantum computation, Phys. Rev. A 79, 022107 (2009).
  * Albash and Lidar [2015] T. Albash and D. A. Lidar, Decoherence in adiabatic quantum computation, Physical Review A 91, 062320 (2015).
  * Zhu  _et al._ [2016] Z. Zhu, A. J. Ochoa, S. Schnabel, F. Hamze, and H. G. Katzgraber, Best-case performance of quantum annealers on native spin-glass benchmarks: How chaos can affect success probabilities, Physical Review A 93, 012317 (2016).
  * Albash  _et al._ [2017] T. Albash, V. Martin-Mayor, and I. Hen, Temperature scaling law for quantum annealing optimizers, Physical Review Letters 119, 110502 (2017).
  * Albash  _et al._ [2019] T. Albash, V. Martin-Mayor, and I. Hen, Analog errors in Ising machines, Quantum Sci. Technol. 4, 02LT03 (2019).
  * Campbell  _et al._ [2017] E. T. Campbell, B. M. Terhal, and C. Vuillot, Roads towards fault-tolerant universal quantum computation, Nature 549, 172 EP (2017).
  * Jordan  _et al._ [2006] S. P. Jordan, E. Farhi, and P. W. Shor, Error-correcting codes for adiabatic quantum computation, Physical Review A 74, 052322 (2006).
  * Lidar [2008] D. A. Lidar, Towards fault tolerant adiabatic quantum computation, Phys. Rev. Lett. 100, 160506 (2008).
  * Young  _et al._ [2013] K. C. Young, M. Sarovar, and R. Blume-Kohout, Error suppression and error correction in adiabatic quantum computation: Techniques and challenges, Phys. Rev. X 3, 041013 (2013).
  * Ganti  _et al._ [2014] A. Ganti, U. Onunkwo, and K. Young, Family of $[[6k,2k,2]]$ codes for practical, scalable adiabatic quantum computation, Phys. Rev. A 89, 042313 (2014).
  * Bookatz  _et al._ [2015] A. D. Bookatz, E. Farhi, and L. Zhou, Error suppression in Hamiltonian-based quantum computation using energy penalties, Physical Review A 92, 022317 (2015).
  * Jiang and Rieffel [2017] Z. Jiang and E. G. Rieffel, Non-commuting two-local hamiltonians for quantum error suppression, Quant. Inf. Proc. 16, 89 (2017).
  * Marvian and Lidar [2017] M. Marvian and D. A. Lidar, Error Suppression for Hamiltonian-Based Quantum Computation Using Subsystem Codes, Physical Review Letters 118, 030504 (2017).
  * Glaetzle  _et al._ [2017] A. W. Glaetzle, R. M. W. van Bijnen, P. Zoller, and W. Lechner, A coherent quantum annealer with Rydberg atoms, Nature Communications 8, 15813 EP (2017).
  * Hamerly  _et al._ [2019] R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli, A. Marandi, T. Onodera, E. Ng, C. Langrock, K. Inaba, T. Honjo, K. Enbutsu, T. Umeki, R. Kasahara, S. Utsunomiya, S. Kako, K.-i. Kawarabayashi, R. L. Byer, M. M. Fejer, H. Mabuchi, D. Englund, E. Rieffel, H. Takesue, and Y. Yamamoto, Experimental investigation of performance differences between coherent Ising machines and a quantum annealer, Science Advances 5, eaau0823 (2019).
  * Tennant  _et al._ [2022] D. M. Tennant, X. Dai, A. J. Martinez, R. Trappen, D. Melanson, M. A. Yurtalan, Y. Tang, S. Bedkihal, R. Yang, S. Novikov, J. A. Grover, S. M. Disseler, J. I. Basham, R. Das, D. K. Kim, A. J. Melville, B. M. Niedzielski, S. J. Weber, J. L. Yoder, A. J. Kerman, E. Mozgunov, D. A. Lidar, and A. Lupascu, Demonstration of long-range correlations via susceptibility measurements in a one-dimensional superconducting Josephson spin chain, npj Quantum Information 8, 85 (2022).
  * Vinci  _et al._ [2016] W. Vinci, T. Albash, and D. A. Lidar, Nested quantum annealing correction, npj Quantum Information 2, 16017 (2016).
  * Pearson  _et al._ [2019] A. Pearson, A. Mishra, I. Hen, and D. A. Lidar, Analog errors in quantum annealing: Doom and hope, npj Quantum Information 5, 1 (2019).
  * Matsuura  _et al._ [2017] S. Matsuura, H. Nishimori, W. Vinci, T. Albash, and D. A. Lidar, Quantum-annealing correction at finite temperature: Ferromagnetic $p$-spin models, Physical Review A 95, 022308 (2017).
  * Trevisan [2004] L. Trevisan, Inapproximability of Combinatorial Optimization Problems, arXiv e-prints (2004), arxiv:cs/0409043 .
  * Arora and Barak [2009] S. Arora and B. Barak,  _Computational Complexity: A Modern Approach_ (Cambridge University Press, 2009).
  * Lucas [2014] A. Lucas, Ising formulations of many NP problems, Front. Phys. 2, 5 (2014).
  * Boothby  _et al._ [2020] K. Boothby, P. Bunyk, J. Raymond, and A. Roy, Next-generation topology of d-wave quantum processors (2020), arXiv:2003.00133 [quant-ph] .
  * Boixo  _et al._ [2013] S. Boixo, T. Albash, F. M. Spedalieri, N. Chancellor, and D. A. Lidar, Experimental signature of programmable quantum annealing, Nat. Commun. 4, 2067 (2013).
  * Bernaschi  _et al._ [2021] M. Bernaschi, M. Bisson, M. Fatica, E. Marinari, V. Martin-Mayor, G. Parisi, and F. Ricci-Tersenghi, How we are leading a 3-XORSAT challenge: From the energy landscape to the algorithm and its efficient implementation on GPUs, Europhysics Letters 133, 60005 (2021).
  * Zurek  _et al._ [2005] W. H. Zurek, U. Dorner, and P. Zoller, Dynamics of a quantum phase transition, Physical Review Letters 95, 105701 (2005).
  * del Campo [2018] A. del Campo, Universal statistics of topological defects formed in a quantum phase transition, Physical Review Letters 121, 200601 (2018).

Methods

Ground state energies
We solved all S28 instances to optimality using Gurobi 10 within feasible
runtime, except for $7$ instances of size $L=15$. For these remaining
instances, Gurobi proved an optimality gap of at most 1.2%, meaning the lowest
energy found is guaranteed to be no greater than $1.2\%$ above the true ground
state. Furthermore, the lowest energies found by PT-ICM were no higher than
those found by Gurobi for all instances. Thus, we assigned the ground state
energies of the remaining instances to the values found by PT-ICM with high
confidence. As we use a median over instances as our summary statistic, our
conclusions are unaffected by the possibility that the ground state energy was
never reached for such few instances.

PT-ICM parameter optimization
We first ran every replica once for $N_{\text{sw}}^{\max}=500,000$ sweeps to
ensure the ground state energy was reached. For the largest size $L=15$, we
determined $N_{\text{sw}}^{(90\%)}$, the number of sweeps that were required
for $90\%$ of the instances to reach their lowest recorded energy, where we
found $N_{\text{sw}}^{(90\%)}\approx 31,000$ for the S28 instances. As this
was significantly less than $N_{\text{sw}}^{\max}$, we considered the ground
states for these instances as “validated” by PT-ICM to calculate median
quantities over the instances. We finally ran PT-ICM 100 times for each
instance, setting $N_{\text{sw}}=N_{\text{sw}}^{(90\%)}$. This yields an
empirical cumulative density function for $p_{E\leq E_{0}+\varepsilon}$ as a
function of the runtime of PT-ICM. The $\operatorname{\text{TT}\varepsilon}$
is then evaluated for each instance by optimizing over the runtime of the PT-
ICM repetition (where $t_{f}$ is now the time needed to reach the target
rather than the annealing time).

The scaling of PT-ICM is ideally evaluated using the parameters that best
optimize the $\operatorname{\text{TT}\varepsilon}$ for each disorder
realization, instance size, and target $\varepsilon$. However, a rigorous
optimization of the number and choice of replica temperatures for all target
$\varepsilon$’s and system sizes is computationally infeasible. To ensure our
results hold for any choice of reasonably optimized temperatures, we repeated
the $\operatorname{\text{TT}\varepsilon}$ evaluation with four temperature
sets summarized in Table 1, which includes both logarithmically-spaced and
feedback-optimized [1] temperatures. The $\operatorname{\text{TT}\varepsilon}$
for a given disorder, instance size, and energy target was chosen from the
best $\operatorname{\text{TT}\varepsilon}$ out of the four temperature sets,
illustrated in Fig. 4. At the optimality gap targets of interest, most
temperature sets’ differences are not appreciable and are unlikely to affect
our conclusions. A more comprehensive but computationally expensive
optimization of PT-ICM would involve a grid search over a range of the
parameters $\beta_{\mathrm{min}}$, $\beta_{\mathrm{max}}$, and
$N_{\mathrm{icm}}$.

| S28
---|---
Set | $N_{T}$ | $\beta_{\mathrm{min}}$ | $\beta_{\max}$ | $N_{\mathrm{icm}}$
1 | 32 | 0.1 | 5.0 | 8
2 | 24 | 0.2 | 10.0 | 6
3 | 32 | 0.2 | 10.0 | 8
4 | 32 | 0.1 | 20.0 | 8
Table 1:  Temperature sets used for PT-ICM for S28 instances, where
$\beta_{\mathrm{min}}$ is the hottest temperature, $\beta_{\max}$ is the
coldest temperature, and $N_{\mathrm{icm}}$ is the number of low temperature
(largest $\beta$) subject to ICM moves. The temperatures are logarithmically
spaced in sets 1 and 4. The temperatures in sets 2 and 3 were feedback-
optimized with initially logarithmically spaced temperatures. Figure 4: Time-
to-epsilon of the different temperature sets used for PT-ICM and their
individual scaling for $\varepsilon=0.75\%$ (left) and $\varepsilon=1.1\%$
(right). The four sets correspond, in order, to the sets shown in Table 1. B5:
Logarithmically-spaced with $\beta_{\mathrm{max}}=5$. F24: Feedback-optimized
with 24 temperatures. F32: Feedback-optimized with 32 temperatures. B20:
Logarithmically-spaced with $\beta_{\mathrm{max}}=20$.

Software Implementation
Our PT-ICM implementation is available as the TAMC software package [2]. Our
implementation of QAC for the D-Wave Advantage graph is available as part of
the PegasusTools Python package [3].

We also implemented and considered the performance of simulated annealing but
do not present its results here as this algorithm was not competitive at large
problem sizes. Our PT-ICM implementation is a general-purpose solver written
in the Rust programming language, targeting CPU execution, and accepts
instances with any connectivity. This is in contrast to previous studies that
have used simulated annealing or the Hamze-de Freitas-Selby algorithm solvers
that are specialized for problems defined on the Chimera graph [4, 5].

Statistical methods
For a given problem size and target optimality gap, we calculate the median
$\operatorname{\text{TT}\varepsilon}$ from the quantum annealing data using a
three-step Bayesian bootstrap procedure over the levels of readout samples,
gauges (spin-reversal transformations at the hardware level), and instances:
(1) the success probability for each gauge is resampled from a beta
distribution for $N_{\text{samp}}$ samples per gauge, (2) the statistical
weight of each gauge is sampled from a Dirichlet distribution of length
$N_{\text{gauge}}$ to take the weighted average success probability for each
instance, and (3) the statistical weight of each instance is sampled from a
Dirichlet distribution of length $N_{\text{inst}}$ to take the weighted
median. We performed our quantum annealing experiments with $N_{\text{samp
}}=1000\mathrm{(QAC)}/3000\mathrm{(U3)}$, $N_{\text{gauge}}=10$, and
$N_{\text{inst}}=125$. We performed $N_{\text{boots}}=200$ bootstrap samples
per size and energy density target pair and found the annealing time and
penalty strength that resulted in an optimal
$\operatorname{\text{TT}\varepsilon}$ for each bootstrap sample. The
distribution of the optimal median $\operatorname{\text{TT}\varepsilon}$
values and the distribution of optimal annealing parameters are the two final
products of this sampling procedure shown in Fig. 2, with $2\sigma$ error
intervals for the optimal $\operatorname{\text{TT}\varepsilon}$.

Next, we formulate a null hypothesis and compute $P$ values as follows. The
null hypothesis is that the optimal annealing time $t_{f}$ is the minimum
accessible $t_{f}^{\min}=0.5\,\mathrm{\mu s}$, i.e., that the true optimal
annealing time is not above $t_{f}^{\min}$. The $P$ value is the empirical
number of bootstrap samples whose optimal $t_{f}$ was $0.5\,\mathrm{\mu s}$,
out of 200 samples. Filled circles in Fig. 2 mean $P<0.05$ for the probability
that $t_{f}=0.5\mu s$ in the bootstrap sample, while open circles mean
$P<0.20$. The filled circles show which points have the highest confidence
that the optimal annealing time is not below $0.5\,\mathrm{\mu s}$.

Non-universality of the observed KZ critical exponent
The Binder cumulant $U$ is well-known to provide a statistical signature of
phase transitions. Under the dynamic finite-size scaling ansatz, $U(L,t_{f})$
is expected to collapse onto a common curve for all system sizes $N$ when
$t_{f}$ is rescaled by $L^{-z-1/\nu}$, where $\nu$ and $z$ are the correlation
length and dynamic critical exponents, respectively. This reflects the KZ
ansatz: the annealing time required for the system to remain adiabatic up to a
correlation length of $L$ scales as

$t_{f}(L)\sim L^{\mu},\quad\mu=z+\frac{1}{\nu}.$ (4)

While the collapse seen in Fig. 3 is visually convincing, the error estimates
for $U$ are, unfortunately, too large to determine the KZ exponent or to
extract $\nu$ and $z$. In addition, we do not observe that the extracted KZ
exponent is universal: the estimate for $\mu$ that best collapses the Binder
cumulant for binomial (as opposed to S28) spin glass instances is $\mu\approx
9\pm 1$ for U3, and is reduced to $\mu\approx 7.5\pm 0.6$ with QAC at
$\lambda=0.1$. Thus, while the scaling ansatz yields a useful complementary
way to quantify the advantage of QAC in terms of sampled quantities in
addition to the $\operatorname{\text{TT}\varepsilon}$ metric, the lack of
universality and the possibility that the spin-glass transition temperature is
zero [6] do not clearly support a universal critical description of the
annealing dynamics, quite unlike the conclusions of Refs. [7, 8].

Time-to-residual energy
Ref. [21] performed energy decay measurements of QA dynamics in 3D spin glass
instances as captured through the residual energy, a dimensionless quantity,

$\rho=\frac{\expectationvalue{H_{z}}-E_{0}}{JN},$ (5)

where $N$ is the number of spins and $J$ is the characteristic energy scale of
the Ising Hamiltonian (which is simply 1 for all of our instance classes).
This motivates an alternative measure for approximate optimization, which we
call _time-to-residual energy_ , or _time-to-rho_ ,

$\mathrm{TTR}(N)=t_{f}\frac{\log(1-0.99)}{\log(1-p_{E\leq E_{0}+\rho JN})}\,,$
(6)

where now $\rho$ sets the target energy difference from the ground state in
units of $JN$, rather than $E_{0}$ as in the case of
$\operatorname{\text{TT}\varepsilon}$. Nevertheless, both $\varepsilon$ and
$\rho$ are targets that grow in proportion to the problem size or number of
variables (in finite dimensions). The TTR analog of Fig. 2 at $\rho=1.1\%$ is
shown in Fig. 5. While the precise scaling exponents $\alpha$ of U3 and QAC
vary slightly, the optimal annealing times are nearly identical. The trends of
the scaling exponents as functions of either $\rho$ or $\varepsilon$ are also
qualitatively similar. In finite-dimensional spin glasses, the variance of
$E_{0}$ across instances appears to scale with $\sqrt{N}$ in the thermodynamic
limit [9]. Hence, $E_{0}/NJ$ converges to a constant, and there is effectively
no distinction between $\varepsilon$ and $\rho$ in the thermodynamic limit.

Figure 5: Time-to-rho performance and scaling of S28 instances, which are
comparable to those using the time-to-epsilon metric.

Classical complexity of approximate optimization
Under the TTR metric, it can be shown that the approximate optimization of
finite-range spin glasses has, in theory, a linear scaling using a simple
divide-and-conquer algorithm. Namely, for a given residual energy target
$\rho$, an algorithm exists whose TTR scales linearly in the system size for a
sufficiently large size, with a large prefactor depending on the residual
energy. For 2D spin glasses on a square lattice, it can be summarized as
follows:

  1. 1.

Partition the size-$N=L^{2}$ instance into $K\times K$ subgraphs $G_{x,y}$,
with $x,y\in\\{1,\ldots,K\\}$, each with _constant_ side length of $L_{0}=L/K$
spins;

  2. 2.

Find the local ground states for each subgraph using an exact or heuristic
solver for each $L_{0}\times L_{0}$ instance. This step has cost $C(L_{0})$,
potentially exponential in $L_{0}^{2}$ on a non-planar lattice;

  3. 3.

Patch together all subgraph ground state configurations as the approximate
ground state for the global Hamiltonian and return this state’s energy as the
approximate ground state energy $E^{*}(N)$. This last step requires $O(N)$
time due to the need to sum the energies of $K^{2}=O(N)$ patches.

The overall complexity is, therefore, $C(L_{0})O(N)$.

This algorithm optimizes _bulk_ energies throughout the volume of the spin-
glass at the expense of large energy violations over regions scaling with the
surface area of the spin-glass, i.e., at the patch boundaries. The boundary
excitations become negligible for sufficiently large sizes compared to the
bulk energies, with the latter eventually reaching the desired residual
energy. More precisely, up to $4L_{0}$ boundary couplings may be violated per
volume of $L_{0}^{2}$, so the residual energy density for this algorithm is
upper-bounded by $4/L_{0}$. Thus, we estimate that the regime of system sizes
where this algorithm applies for the $\rho$ targets we examine, e.g., 1%,
would require a reliable, efficient solver for instances with at least a _sub-
problem_ side length of $L_{0}\approx 400$, resulting in a prefactor of
$2^{400}\sim 10^{120}$, which is entirely impractical. Nevertheless, this
algorithm could be a starting point for parallel and quantum-classical hybrid
algorithms for massive, finite-dimensional problems. Furthermore, such an
algorithm could reach a similar target $\varepsilon$ in linear time with
increasing probability as the system size increases due to the decreasing
variance of $E_{0}/(JN)$.

Figure 6:  Time-to-epsilon scaling of spin glass instances with R6 disorder
for (a) $\varepsilon=0.75\%$ and (b) $\varepsilon=1.10\%$. To expedite finding
the optimal $\operatorname{\text{TT}\varepsilon}$, the annealing times were
limited to the range $t_{f}\in[0.5,9]\,\mathrm{\mu s}$ and QAC was limited to
$J_{p}=0.2$ only.

Alternative spin-glass disorder cases
We performed a similar study of instances with binomial disorder $J=\pm 1$.
UnsurprIsingly, there was no substantial difference in scaling between U3 and
QAC, and additionally, we were unable to validate the optimal annealing time
for QAC at $\varepsilon=1\%$. Thus, we considered the binomial disorder
instances unsuitably easy for our purposes.

In an attempt to estimate the influence of precision and ground state
degeneracy in the S28 disorder, we also studied range 6 (R6) instances for
which the couplings were randomly drawn from $J\in\\{\pm 1/6,\pm 2/6,\pm
3/6,\pm 4/6,\pm 5/6,\pm 1\\}$. The main contrast between R6 and S28 is that
the minimum non-zero local field a qubit may experience is much smaller in S28
disorder. Furthermore, in the absence of the Sidon property, the R6 disorder
is more susceptible to ground-state degeneracies due to floppy qubits.
However, it is less likely to occur than in a smaller range disorder case such
as range 3.

The $\operatorname{\text{TT}\varepsilon}$ scaling for the R6 case is shown in
Fig. 6. For a smaller $\varepsilon=0.75\%$, the scaling of QAC is better with
R6 disorder than S28 disorder (Fig. 2), going from 2.26 to 2.02, then
equalizes as $\varepsilon$ increases. Perhaps surprIsingly, the opposite is
true of U3, which scales worse on R6 disorder at smaller epsilon (2.92 to
3.10) before scaling better than S28 disorder above $\varepsilon\approx 1.1\%$
(2.02 to 1.82). That is, while S28 is overall more challenging than R6 under
QAC, low energy sampling of R6 is more challenging for unprotected QA than
S28. This is likely caused by the worse relative precision of small coupling
strengths in unprotected QA, which would not affect the S28 disorder.

## References

  * Katzgraber  _et al._ [2006] H. G. Katzgraber, S. Trebst, D. A. Huse, and M. Troyer, Feedback-optimized parallel tempering Monte Carlo, Journal of Statistical Mechanics: Theory and Experiment 2006, P03018 (2006).
  * H. Munoz Bauza [2023a] H. Munoz Bauza, TAMC software package (2023a).
  * H. Munoz Bauza [2023b] H. Munoz Bauza, PegasusTools Python package (2023b).
  * Hamze and de Freitas [2004] F. Hamze and N. de Freitas, From fields to trees, in  _Proceedings of the 20th Conference on Uncertainty in Artificial Intelligence_ , UAI ’04 (AUAI Press, Arlington, Virginia, USA, 2004) pp. 243–250.
  * Selby [2014] A. Selby, Efficient subgraph-based sampling of Ising-type models with frustration, arXiv:1409.3934 (2014).
  * Katzgraber  _et al._ [2014] H. G. Katzgraber, F. Hamze, and R. S. Andrist, Glassy Chimeras could be blind to quantum speedup: Designing better benchmarks for quantum annealing machines, Physical Review X 4, 021008 (2014).
  * King  _et al._ [2022] A. D. King, S. Suzuki, J. Raymond, A. Zucca, T. Lanting, F. Altomare, A. J. Berkley, S. Ejtemaee, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, J. D. Whittaker, J. Yao, R. Harris, D. A. Lidar, H. Nishimori, and M. H. Amin, Coherent quantum annealing in a programmable 2,000 qubit Ising chain, Nature Physics 18, 1324 (2022).
  * King  _et al._ [2023] A. D. King, J. Raymond, T. Lanting, R. Harris, A. Zucca, F. Altomare, A. J. Berkley, K. Boothby, S. Ejtemaee, C. Enderud, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, R. Molavi, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, N. Tsai, M. Volkmann, J. D. Whittaker, J. Yao, A. W. Sandvik, and M. H. Amin, Quantum critical dynamics in a 5,000-qubit programmable spin glass, Nature 617, 61 (2023).
  * Bouchaud  _et al._ [2003] J.-P. Bouchaud, F. Krzakala, and O. C. Martin, Energy exponents and corrections to scaling in Ising spin glasses, Physical Review B 68, 224404 (2003).