# Dependence of Lightcurves on Phase Angle and Asteroid Shape

Xiao-Ping Lu Faculty of Information Technology, Macau University of Science
and Technology, Taipa, Macao, China  State Key Laboratory of Lunar and
Planetary Sciences(Macau University of Science and Technology)  David Jewitt
Department of Earth, Planetary and Space Sciences, UCLA, 595 Charles Young
Drive East, Los Angeles, CA 90095-1567, USA  Department of Physics and
Astronomy, UCLA, 430 Portola Plaza, Box 951547, Los Angeles, CA 90095-1547,
USA

(Accepted in AJ, 2019 October 8)

###### Abstract

We investigate the phase angle dependence of asteroid lightcurves using
numerical scattering models applied to simple body shapes. For simplicity, the
Kaasalainen scattering parameters are obtained from the corresponding Hapke
scattering laws for C-type, S-type, and V-type asteroids. The spectral types
differ substantially in the role of multiple scattering (which is largely a
function of their geometric albedos) but we find that the differences on the
light curve vs. phase relations are modest. By using a Kaasalainen scattering
law, the amplitudes and axis ratios with respect to different phase angles
from $0^{\circ}$ to $140^{\circ}$ are plotted for these types asteroids based
on a biaxial ellipsoid shape model. Additionally, we examine the relationship
between amplitude and axis ratio for a contact binary represented by identical
biaxial ellipsoids, including the effects of shadowing of one component by the
other. We compare the models with published high phase angle observations, and
with interstellar object 1I/‘Oumuamua ($\Delta m=2.5$ magnitude at phase angle
$\alpha=23^{\circ}$), finding axis ratio $5.2:1$ if represented as a single
ellipsoid and $3.5:1$ for each component if represented as a nose-to-nose
contact binary. While a detailed fit is not attempted, the comparison shows
that the single ellipsoid model is better.

techniques: photometric – methods: numerical – minor planets, asteroids:
general

##  1 Introduction

A huge number of photometric observations of asteroids are obtained everyday,
greatly enriching our knowledge of this population. The most basic parameters
to be extracted from photometric data include estimates of the asteroid size
(if the albedo is known or can be assumed) and the shape. Of course, while the
true shapes of asteroids are irregular, it is convenient to represent them by
geometric forms, typically as triaxial ellipsoids characterized by their axis
ratios. The axis ratios can then be inferred from observed brightness
variations as the body rotates.

However, in addition to the intrinsic shape, it is well known that the light
curve shape and amplitude also depend on the angle between the spin vector and
the line of sight, and on the phase (Sun-object-observer) angle, $\alpha$.
Generally the amplitude grows larger (because of scattering and self-shadowing
effects), for a given axis ratio, as the phase angle increases. Zappalà et al.
(1990) first analyzed the relationship between light curve amplitude and solar
phase angle among asteroids. They found that the amplitude is linearly
proportional to phase angle for $\alpha\leq 50^{\circ}$. The phase angle
dependence of the light curve has also been investigated by Lacerda & Jewitt
(2007) but has otherwise not received much quantitative attention.

In this paper we investigate the relationship between amplitude, axis ratio,
spectral type and phase angle based on biaxial ellipsoid models, in order to
provide a quick and practical tool for observers. The numerical method of
calculating the amplitude is based on a biaxial ellipsoid shape with the Hapke
scattering law and Kaasalainen law, as presented in Section 2. The
relationship between amplitude and axis ratio is discussed in Section 3,
together with a comparison of phase effects in different asteroid spectral
types. We briefly consider applications to data in Section 4.

##  2 Scattering Law

###  2.1 Simulated Brightness

To invert photometric observations (Kaasalainen & Torppa, 2001; Muinonen &
Lumme, 2015), the simulated brightness is calculated from the surface integral
based on different shape models, scattering laws, and viewing geometries,
using

$L=\iint_{C+}S(\alpha,\mu,\mu_{0})\,d\sigma.$ (1)

Here, $S(\alpha,\mu,\mu_{0})$ is the scattering function and $C+$ represents
the part of surface both illuminated and visible for different shape models,
commonly considered as triaxial ellipsoids (Cellino et al., 2015), convex
polyhedron shapes (Kaasalainen et al., 2012), and the intermediate Cellinoid
shapes (Lu et al., 2016, 2017). $\mu,\mu_{0}$ in Equation (1) denote the
projections of unit vectors of viewing and illuminating directions on the
surface normal direction.

###  2.2 Scattering Functions

Hapke (2012) introduced a semi-physical model for illustrating the light
reflection behaviour of a planetary surface. Here we exploit the five-
parameter version for the scattering law, $S_{HM}(\alpha,\mu,\mu_{0})$,
written as

$S_{HM}(\alpha,\mu,\mu_{0})=\frac{\omega}{4}\frac{\mu_{0}}{\mu_{0}+\mu}\left\\{\left[1+B(\alpha)\right]p(\alpha)+H(\mu_{0})H(\mu)-1\right\\}S(i,e,\alpha,\bar{\theta}),$
(2)

where the five parameters are the single scattering albedo ($\omega$),
roughness ($\bar{\theta}$), asymmetry factor ($g$), opposition surge amplitude
($B_{so}$), and opposition surge width ($h_{s}$), respectively. The related
formulae, such as shadow hiding opposition surge function $B(\alpha)$, average
particle single-scattering phase function $p(\alpha)$, and Ambartsumian-
Chandrasekhar function for multiple scattering $H(x)$ can be found in Hapke
(2012) (page 323).

For simplicity, some empirical scattering functions are frequently used in
numerical simulations. For example, Kaasalainen et al. (2001) used a three-
parameter empirical phase function, which depends only on phase angle, with an
empirical disk function (which depends solely on incidence and emission
angles) to describe the surface scattering function,
$S_{MK}(\alpha,\mu,\mu_{0})$. The phase function is given by

$f(\alpha)=A_{0}\exp\left(-\frac{\alpha}{D}\right)+K\alpha+1,$ (3)

where the opposition effect is described by a simple exponential term with
height ($A_{0}$) and width ($D$) followed by a linear term with slope, $K$, to
describe the phase curve behavior outside of opposition. The phase function is
a combination of a single (Lommel-Seeliger function) and multiple (Lambert
function) term, with a relative weighting factor, $C$, such that,

$S_{MK}(\alpha,\mu,\mu_{0})=f(\alpha)\left(\frac{\mu\mu_{0}}{\mu+\mu_{0}}+C\mu\mu_{0}\right).$
(4)

###  2.3 Asteroid Taxonomy

Hapke photometric parameters representing asteroids of different spectral
types are listed in Table 1. Li et al. (2015) listed the corresponding five
parameters of the Hapke model for S-type, C-type, and V-type asteroids. The
C-type asteroids have an average albedo $\omega=0.037$, indicating a dark
surface, while V-type asteroids have a brighter average albedo $\omega=0.51$.
The S-type asteroids are intermediate, with the average albedo $\omega=0.23$.

The Hapke model is difficult to invert numerically and so, instead, we fit the
related parameters of the Kaasalainen model (MK Parameters) for different
asteroid spectral types. First, the simulated light curve ($LC_{HM}$) is
generated based on one specified ellipsoid shape ($C$) with the supposed
viewing($E$) and illuminating($E_{0}$) geometries, as well as the Hapke
scattering parameters ($P_{HM}$) in Equation (2) for some taxonomy as listed
in Table 1. Each point in the light curve ($LC_{HM}$) corresponding to one
rotational phase angle can be calculated by the surface integration using
Equation (1) as follows,

$LC^{(i)}_{HM}(C,E,E_{0},P_{HM})=\iint_{C+}S_{HM}(\alpha,\mu,\mu_{0})\,d\sigma.$
(5)

Then we search for the best-fit parameters of the Kaasalainen model by
minimizing,

$\chi^{2}=\||LC_{HM}-LC_{MK}\||_{2}^{2}$ (6)

where $LC_{MK}$ denotes the simulated light curve generated based on MK
scattering parameters in Equation (4). Each point in $LC_{MK}$ is calculated
under the same conditions as $LC_{HM}$, except for the scattering model,

$LC^{(i)}_{MK}(C,E,E_{0},P_{MK})=\iint_{C+}S_{MK}(\alpha,\mu,\mu_{0})\,d\sigma.$
(7)

It should be noted that as the simulated lightcurves here are not relative
(i.e. Normalization), the general form of the phase function is exploited as
follows,

$f(\alpha)=A_{0}\exp\left(-\frac{\alpha}{D}\right)+K\alpha+B.$ (8)

with the additional back-ground intensity $B$ in Kaasalainen model, as
presented in Kaasalainen et al. (2003).

Finally, we test more than 1000 cases covering different shape models and
observing geometries to obtain the distributions of MK parameters for each
spectral type. The obtained MK parameters are listed in Table 1, together with
the Hapke parameters for comparison.

Table 1 shows that as the albedo increases, the weight factor $C$ also
increases, meaning that multiple scattering effects increase as the albedo
increases (consistent with Hapke (2012)). Furthermore, the parameter $C$ in
the Kaasalainen scattering law is evidently related to the asteroid taxonomy.
It can be applied to roughly classify the asteroids based on the derived MK
parameters from photometric observation. In the numerical simulation, the
appropriate $C$ value can improve the accuracy of other parameters, such as
the pole direction and shape, if the asteroid types are known from other
measurements, such as infrared observations (Wright et al., 2010).

Furthermore, multiple scattering is intrinsic to the asteroid surface and,
unlike the phase function, is not cancelled when calculating the amplitude.
Therefore, for simplicity in the following section the relationship between
amplitude and axis ratio for different taxonomical asteroids will adopt the
derived $C$ in Table 1.

##  3 Amplitude and Axis Ratio

The observed brightness of an asteroid is related to the viewing and
illumination geometries, as well as to its physical parameters, including pole
orientation and overall shape. Here we suppose that the shape is a biaxial
ellipsoid ($a>b=c$) in rotation about the $c$-axis, and the viewing and
illuminating directions are located on the equatorial plane of the asteroid.
With this simplification, the simulated amplitude will provide an estimate of
the upper limit to the amplitudes in different circumstances.

The lightcurve amplitude $\Delta m$ is calculated from

$\Delta m=2.5\log\left(\frac{L_{\max}}{L_{\min}}\right)$ (9)

where $L_{\min}$ and $L_{\max}$ are the minimum and maximum values of the
calculated lightcurve, respectively.

The Kaasalainen scattering law in Equation (4) is applied to calculate the
simulated brightness in Equation (1). The primary contributions to the
amplitude, $\Delta m$, are the disk function with the weight factor $C$ and
the shape.

We computed the relationships between lightcurve amplitude, $\Delta m$, and
axis ratio, $(a/b)$, with respect to phase angle for C-type, S-type, and
V-type asteroids. The results for C-types are shown in Figure 2 while those
for S-types are in Figure 4. The results for V-types are very similar and so
are not separately plotted here. Weidenschilling (1980) noted that a rotating
fluid (strengthless) body cannot adopt an equilibrium shape if $a/b>$ 2.3.
Indeed, almost all well-observed asteroids have $a/b$ smaller than this value
but, at small sizes where strength effects are important, larger ratios are
possible. With this in mind, we consider a range of axis ratios from $a/b$ =
2:1 to $a/b$ = 7:1. We extend the models in Figures 2 and 4 only to $\alpha$ =
140° because, at larger phase angles, the computed lightcurves become very
sensitive to the adopted shape and scattering function and, in any case,
essentially no observations exist at such large angles.

For small values of $a/b$, Figures 2 and 4 show that $\Delta m\propto\alpha$
for $\alpha\gtrsim 20$°, broadly confirming the result of Zappalà et al.
(1990). At phase angles $\alpha\lesssim 20^{\circ}$, the amplitude is a
slightly non-linear function of $\alpha$, as a result of the opposition
effect. As the phase angle approaches $0^{\circ}$, the observed brightness
will increase exponentially. Comparison of Figures 2 and 4 shows that the
phase angle dependence of the amplitudes, for a given shape, is only weakly
dependent on the asteroid spectral type. However, apparently the amplitude
slightly increases as the albedo increases, especially as observed at the
small phase angles.

##  4 Applications

An obvious conclusion from Figures 2 and 4 is that high phase angle
observations will preferentially produce large amplitude lightcurves.
Measurements of the shapes and shape distributions of asteroids, if based on
observations taken at non-zero phase angles, should be interpreted with
appropriate caution. Likewise, rotational lightcurves accumulated over a wide
range of phase angles can be compromised by the phase angle dependence, as is
apparently the case in recent spacecraft observations of Kuiper belt objects
2012 HE85 and 2011 HK103 (Verbiscer et al., 2019).

A less obvious conclusion is that the exaggerated lightcurves will lead to a
size-dependent bias in asteroid surveys undertaken at large phase angles (e.g.
the small elongation/large phase angle “sweet spot” survey strategy originally
proposed for Panstarrs; Wainscoat et al. 2015). Most asteroid sky surveys are
magnitude-limited and objects whose mean magnitude is close to the survey
detection limit will oscillate in and out of detection, according to the
instantaneous rotational phase. The result is a bias against the detection of
faint asteroids, because they will be missing in many of the repeat visits
needed by automatic detection routines for their identification. For example,
at phase angle $\alpha$ = 100°, a body with unremarkable axis ratio $a/b$ =
2:1 would show a lightcurve amplitude $\Delta m\sim$ 1.8 magnitudes (a factor
of 5). Objects having a mean brightness within $\Delta m$ magnitudes of the
survey limit would dip into and out of detection according to the rotational
phase. Bright asteroids (i.e. those with mean magnitude brighter than the
survey limit by an amount $>\Delta m$) can remain above the detection
threshold regardless of rotational phase. Since brightness and asteroid size
are (statistically) related, the under-counting of faint asteroids will lead
to an apparent flattening of the size distribution, and an under-estimate of
the flux of Earth-approaching bodies in asteroid surveys when conducted away
from opposition.

The models show that the angle dependence of the lightcurve amplitude is
slightly influenced by multiple scattering and is therefore related to the
albedo and so to the spectral type. In principle, this fact could be used to
spectrally classify asteroids even using single filter measurements taken over
a range of geometries. Such measurements might be produced, for example, by
wide field sky surveys taken for other scientific purposes. However, the
spectral type dependence is very weak (compare Figures 2 and 4) so that high
quality observations would be needed for this single filter spectral-typing to
be possible.

To date, few high phase angle lightcurves appear in the literature. To give
one example, Denk et al. (2018) presented observations of the irregular
satellites of Saturn, including some high phase angle lightcurves. The two
objects with the largest lightcurve amplitudes are plotted in Figure 6.
Lightcurves of 15 km scale Kiviuq were determined at phase angle
$\alpha=108^{\circ}$ with amplitude $\Delta m=2.5$ magnitude, and 40 km scale
Siarnaq was observed at $\alpha=121^{\circ}$ with $\Delta m=2.0$ magnitude,
respectively. These very large amplitudes, if interpreted geometrically, would
indicate axis ratios 100.4Δm = 10:1 and 6:1, respectively. Inspection of
Figure 6 shows that inclusion of the phase angle dependence of the scattering
gives a much more modest axis ratios from $\sim$2.0 to 2.6.

As another example, the so-called interstellar asteroid, 1I/(2017 U1)
‘Oumuamua (hereafter “U1”), displayed an abnormally large rotational
amplitude, $\Delta m\sim 2.5$ magnitude that was interpreted as axis ratio
$a/b\geq$ 10 by Meech et al. (2017). However, U1 was not observed at
$\alpha=0$° and so the phase angle vs amplitude effect must be taken into
account. Specifically, observations of U1 were possible over only a very
narrow range of dates, with most conducted between October 25 and October 30,
2017 when the solar phase angle was $\alpha\sim 23^{\circ}$. The albedo of
this object could not be determined (Trilling et al., 2018), but a comet-like
geometric albedo $p=0.04$ is widely assumed (Meech et al., 2017; Fitzsimmons
et al., 2018). With $\alpha=23^{\circ}$, $\Delta m$ = 2.5 magnitudes and $p$ =
0.04 (corresponding to a C-type surface) we can estimate from Figure 2 an axis
ratio $a/b=5$. This result is similar to the estimate by Bannister et al.
(2017), who gave an estimate for the axis ratio about $5.3:1$. This is still a
much larger axis ratio than found among the main-belt asteroids, most of which
have $a/b<$ 2, but less extreme than obtained when neglecting the phase angle
effect.

As an experiment, we also simulated the lightcurve for an assumed contact
binary shape in which the two components are identical biaxial ellipsoids
(c.f. Lacerda and Jewitt (2007), who used Jacobi ellipsoids and different
scattering functions). The relationship between amplitude and axis ratio is
shown in Figure 7, where we have assumed a C-type surface with albedo
$p=0.037$. In this case, we find that the lightcurve of U1 can be approximated
by nose-to-nose biaxial ellipsoids each with axis ratio $a/b=3.5$, as marked
in the figure.

Figure 8 shows simulated lightcurves from the two different shape models,
overplotting the lightcurve data from Drahus et al. (2018). In the top panel
we plot the single biaxial ellipsoid model with axis ratio $a/b=5.2$ while the
lower panel shows a contact binary model having equal components each with
$a/b=3.5$. Both shape models can crudely reproduce the lightcurve amplitudes.
We have made no attempt to match the detailed shape of the lightcurve because,
in such a small body, the effects of irregularity are likely to be
significant. However, one feature of the lightcurve that is less sensitive to
irregularity or other details of the shape is the flat minimum present in the
contact binaries caused by mutual shadowing (Figure 8 and Lacerda and Jewitt
2007). The lightcurve of U1 does not show a flat minimum. For this reason, we
prefer the single ellipsoid model, and conclude that U1 has an axis ratio
closer to $a/b\sim$ 5 than to $a/b\sim$ 10, in agreement with Drahus et al.
(2018).

##  5 Summary

  * •

We present models of the phase angle dependence of the lightcurve amplitude
for asteroids having different axis ratios, spectral types, and ellipsoid and
contact binary body-shapes. Our results (Figures 2 and 4) provide a quick
constraint on the shape for observations taken at non-zero phase angles.

  * •

We find that, to first order, the amplitude vs. phase relation is independent
of the assumed asteroid albedo and spectral type, and depends only on the body
shape.

  * •

Using the observations of 1I/2017 U1 ‘Oumuamua tabulated by Drahus et al.
(2018), we infer axis ratio $5.2:1$ for an ellipsoid shape model. We can also
match the amplitude with a contact binary model in which each component has
axis ratio $3.5:1$. However, the lightcurve of U1 lacks the flat-bottomed
minima characteristic of contact binaries and so we prefer the former
representation.

We thank the anonymous referees for their constructive comments. XL is funded
by The Science and Technology Development Fund, Macau SAR (File no.
0018/2018/A).

## References

  * Bannister et al. (2017) Bannister, M. T., Schwamb, M. E., Fraser, W. C., et al. 2017, ApJL, 851, L38
  * Cellino et al. (2015) Cellino, A., Muinonen, K., Hestroffer, D., & Carbognani, A. 2015, Planet. Space Sci., 118, 221
  * Denk et al. (2018) Denk T., Mottola S., Tosi F., Bottke W.F., and Hailton D. P. 2018, The irregular satellites of Saturn. In Enceladus and the Icy Moons of Saturn(P. M. Schenk et al., eds.) (Univ. of Arizonna, Tucson), 409–434
  * Drahus et al. (2018) Drahus, M., Guzik, P., Waniak, W., et al. 2018, Nature Astronomy, 2, 407
  * Ďurech et al. (2010) Ďurech, J., Sidorin, V., & Kaasalainen, M. 2010, A&A, 513, A46
  * Feng & Jones (2018) Feng, F., & Jones, H. R. A. 2018, ApJL, 852, L27
  * Fitzsimmons et al. (2018) Fitzsimmons, A., Snodgrass, C., Rozitis, B., et al. 2018, Nature Astronomy, 2, 133
  * Hapke (2012) Hapke, B. 2012, Theory of reflectance and emittance spectroscopy, 2nd edn. (Wiley Online Library)
  * Jewitt et al. (2017) Jewitt, D., Luu, J., Rajagopal, J., et al. 2017, ApJ, 850
  * Kaasalainen et al. (2012) Kaasalainen, M., Lu, X.-P., & Vänttinen, A. 2012, A&A, 539, A96
  * Kaasalainen & Torppa (2001) Kaasalainen, M., & Torppa, J. 2001, Icarus, 153, 24
  * Kaasalainen et al. (2001) Kaasalainen, M., Torppa, J., & Muinonen, K. 2001, Icarus, 153, 37
  * Kaasalainen et al. (2003) Kaasalainen, S., Piironen, J., Kaasalainen, M., et al. 2003, Icarus, 161, 34
  * Lacerda & Jewitt (2007) Lacerda, P., & Jewitt, D. C. 2007, AJ, 133, 1393
  * Li et al. (2015) Li, J.-Y., Helfenstein, P., Buratti, B. J., Takir, D., & Clark, B. E. 2015, Asteroid photometry. In Asteroids IV (P. Michel et al., eds.) (Univ. of Arizona, Tucson), 129–150
  * Lu et al. (2016) Lu, X.-P., Cellino, A., Hestroffer, D., & Ip, W.-H. 2016, Icarus, 267, 24
  * Lu et al. (2017) Lu, X.-P., Ip, W.-H., Huang, X.-J., & Zhao, H.-B. 2017, Planet. Sp. Sci., 135, 74
  * Meech et al. (2017) Meech, K. J., Weryk, R., Micheli, M., et al. 2017, Nature, 552, 378
  * Moro-Martín (2018) Moro-Martín, A. 2018, ApJ, 866, 131
  * Muinonen & Lumme (2015) Muinonen, K., & Lumme, K. 2015, A&A, 584, A23
  * Trilling et al. (2018) Trilling, D. E., Mommert, M., Hora, J. L., et al. 2018, AJ, 156, 261
  * Verbiscer et al. (2019) Verbiscer, A. J., Porter, S., Benecchi, S. D., et al. 2019, AJ, 158, 123
  * Wainscoat et al. (2015) Wainscoat, R. J., Chambers, K., Lilly, E., et al. 2015, AAS/Division for Planetary Sciences Meeting Abstracts #47, 301.07
  * Weidenschilling (1980) Weidenschilling, S. J. 1980, Icarus, 44, 807
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Ye et al. (2017) Ye, Q.-Z., Zhang, Q., Kelley, M. S. P., & Brown, P. G. 2017, ApJL, 851, L5
  * Zappalà et al. (1990) Zappalà, V., Cellino, A., Barucci, M. A., M., F., & D., L. F. 1990, A&A, 231, 548

Table 1: Derived MK Parameters from Hapke Model for different asteroids Parameters | C-Type | S-Type | V-Type
---|---|---|---
Hapke Model
$\omega$ | 0.037 | 0.23 | 0.51
$g$ | 1.03 | 1.6 | 1
$B_{so}$ | 0.025 | 0.08 | 0.098
$h_{s}$ | -0.47 | -0.27 | -0.26
$\bar{\theta}$ | 20 | 20 | 32
Kaasalainen Model
$A_{0}$ | $0.009\sim 0.01$ | $0.039\sim 0.041$ | $0.046\sim 0.048$
$D$ | $3.3\sim 3.5$ | $8.4\sim 8.6$ | $9.2\sim 9.45$
$K$ |  $-2.4\sim-2.2$ $\times 10^{-4}$ |  $-4.8\sim-4.6$ $\times 10^{-4}$ |  $-7.7\sim-7.3$ $\times 10^{-4}$
$B$ | $0.0142\sim 0.0148$ | $0.049\sim 0.051$ | $0.081\sim 0.087$
$C$ | $0.08\sim 0.11$ | $0.11\sim 0.15$ | $0.28\sim 0.34$

Figure 1:  Figure 2: Lightcurve amplitude, $\Delta m$ vs. phase angle
$0^{\circ}\leq\alpha\leq 140^{\circ}$ for biaxial ellipsoid ($a>b=c$) models
with $a/b$ =2,3…7, and with C-type scattering, assumed average albedo
$p=0.037$.

Figure 3:  Figure 4: Lightcurve amplitude, $\Delta m$ vs. phase angle
$0^{\circ}\leq\alpha\leq 140^{\circ}$ for biaxial ellipsoid ($a>b=c$) models
with $a/b$ =2,3…7, and with S-type scattering, assumed average albedo
$p=0.23$.

Figure 5:  Figure 6: Lightcurve amplitude, $\Delta m$ vs. phase angle for
C-type scattering, with the notations of 2017U1 and Saturn satellites, Kiviuq
and Siarnaq, corresponding the observed phase angles.

Figure 7: Same as Figure 6 but for a contact binary model with C-type
scattering function.

Figure 8: Simulated lightcurves (green line) for U1 at $\alpha$ = 23° based on
single and binary shape models, respectively, with the comparison to the
observations (denoted by red and blue stars) by Drahus et al. (2018).


# Catalysis in Charging Quantum Batteries

B. Ahmadi<EMAIL_ADDRESS>International Centre for Theory of Quantum
Technologies, University of Gdansk, Jana Bażyńskiego 1A, 80-309 Gdansk, Poland
P. Mazurek  International Centre for Theory of Quantum Technologies,
University of Gdansk, Jana Bażyńskiego 1A, 80-309 Gdansk, Poland  R. R.
Rodríguez  International Centre for Theory of Quantum Technologies, University
of Gdansk, Jana Bażyńskiego 1A, 80-309 Gdansk, Poland  S. Barzanjeh
Department of Physics and Astronomy, University of Calgary, Calgary, AB T2N
1N4 Canada  R. Alicki  International Centre for Theory of Quantum
Technologies, University of Gdansk, Jana Bażyńskiego 1A, 80-309 Gdansk, Poland
P. Horodecki  International Centre for Theory of Quantum Technologies,
University of Gdansk, Jana Bażyńskiego 1A, 80-309 Gdansk, Poland

###### Abstract

We demonstrate that energy transfer limitations in the context of charging
quantum batteries can be significantly enhanced in the presence of catalyst
qubits. A catalyst qubit is added to the process of charging to enhance the
amount of energy transferred to the quantum battery (arbitrary amount of
energy transfer!), without consuming almost no energy by the catalyst qubit
throughout the whole process. The use of a catalyst qubit also renders the
setting autonomous, such that no external measurement of the coupling strength
are needed to operate within the optimal energy transfer mode. The charger is
driven by an external laser field. Both the charger and the battery are
modeled by quantum harmonic oscillators.

##  I Introduction

With the advent of quantum technology a growing curiosity, regarding how the
current technologies can be enhanced under the principles of quantum
mechanics, began. Different topics have been explored by scientists, such as
quantum thermodynamics, quantum cryptography, quantum metrology and quantum
computing. The quantum thermodynamics theory addresses the concepts of energy,
work, heat, temperature, and entropy production in quantum systems Spohn
(1978); Alicki (1979); Parrondo et al. (2015); Brandao et al. (2015); Ahmadi
et al. (2021, 2019). Recently the study of energy transfer has brought forward
a diverse and rich set of phenomena, in the context of quantum batteries, to
be both theoretically and experimentally investigated Alicki and Fannes
(2013); Campaioli et al. (2017, 2018); Barra (2019); Farina et al. (2019);
Julià-Farré et al. (2020). As in classical batteries the aim is to boost the
energy transfer, from the battery to the charger, in order to be later
extracted in the form of work. Many authors have proposed and investigated the
superiority and possible benefits of quantum effects, such as quantum
correlations and quantum coherence of the initial state, in charging quantum
batteries Ferraro et al. (2018); Andolina et al. (2019a, b).
Even though this field is promising and filled with a plethora of interesting
applications, there exist some issues or experimental limitations to be
addressed. One of them is creating the initial quantum coherence in the state
of the charger for boosting the energy transfer. This, of course, requires
specific preparation which is not experimentally favorable thus, in general,
it would be preferable to start from the ground state. Quantum correlations
may also be troublesome since uncorrelating the battery from the charger can
significantly disturb the state of the battery and cost some energy as well.
Another issue to be noticed is that it is energetically and experimentally
more beneficial to use quantum harmonic oscillators (QHO), especially as a
battery for storing energy, rather than two-level systems (TLS). The reason is
the ability to store an arbitrary amount of energy (up to natural technical
limitations of the particular experimental setup) in one QHO while the maximum
energy stored by a TLS is restricted to one quanta. To overcome this
limitation, we could use many TLSs. However, this is not experimentally
desirable because of the noise produced when adding more TLSs.
Catalyst qubit was first proposed by scientists in the resource theory for
quantum state transformation purposes such that it allowed for some state
transformations that were impossible under some specific operations Jonathan
and Plenio (1999); Eisert and Wilkens (2000). No operational use of catalyst
qubit has been proposed in the context of energy transfer in quantum
thermodynamic systems. Can they also be used to facilitate more energy
transfer between two quantum thermodynamic systems?

Figure 1: (Color online) Top: efficient charging of a quantum system (battery)
cannot be achieved by external energy driving in the presence of noise.
Bottom: addition of a qubit ancilla makes the energy transfer effective, while
keeping the state of the ancilla almost unaffected.

In this work, we demonstrate that the obstacles mentioned above may be
overcome by the use of catalyst qubit in the process of charging. More
precisely, we will show that the mere presence of a qubit in ground state can
allow charger-battery systems to realize arbitrary energy transfers that would
otherwise be impossible. We first examine the basic model of quantum charger-
battery system. As stated above, TLSs are not appropriate for storing an
arbitrary amount of energy thus, to this end, we use QHOs as both quantum
charger and quantum battery. An external laser field is shined on the charger
and both the charger and the battery are assumed to be initially in the ground
state. We will then examine the intuitive way of boosting energy transfer in
the basic model which is based on the tuning of the external laser field with
the system of the charger-battery. Then we will introduce a catalytic approach
to the task of charging quantum batteries, with an ancillary system added to
the system of the charger-battery. We call it a catalyst qubit because the
qubit does not consume (waste) almost any energy during the charging process
and it remains close to its initial state (ground state) after the process.
More importantly the catalyst does not get significantly entangled with
neither the charger nor the battery. The catalytic approach has the advantage
of being experimentally friendly because one does not need to know about the
interaction strength between the charger and the battery.

##  II Energy transfer from a charger to a battery: The Basic Model

Figure 2: (Color online) A QHO acting as a quantum charger (middle), powered
by the classical laser field (left) and interacting with another QHO acting as
a battery (right). Energy dissipation, indicated by blue curved arrows, occurs
through the charger.

In this section we present the basic model of charger-battery system and in
the next sections we will present how the energy transfer can be greatly
booted in this basic model. Consider a QHO, as a charger, which interacts
(charges) another QHO (as the battery) (see Fig. (2)). A classical laser field
with frequency $\omega$ feeds energy into the charger. Thus the Hamiltonian is
written as ($\hbar=1$)

$\displaystyle H$
$\displaystyle=\omega_{a}a^{\dagger}a+\omega_{b}b^{\dagger}b+g(ab^{\dagger}+ba^{\dagger})$
$\displaystyle+F\left(e^{-i\omega_{f}t}a+e^{i\omega_{f}t}a^{\dagger}\right),$
(1)

where $a$ ($a^{\dagger}$) and $b$ ($b^{\dagger}$) are annihilation (creation)
operators of the charger $A$ and the battery $B$, respectively. $g$ is the
coupling constant, $\omega_{f}$ is the laser field frequency, and $F$ the
amplitude of the laser field. The interaction between the laser field and the
charger has been computed in dipole approximation and rotating wave
approximation (RWA) Loudon and Petruccione (2000). In quantum optics
applications the frequencies are high in comparison to decay rates therefore
RWA (or secular approximation at the level of Master equation) is valid. Both
the charger and the battery are prepared initially in ground states. It is
assumed that the dissipation of energy, into the environment $E$, occurs
through the charger. Therefore the evolution of the state of the whole system
$\rho_{AB}$, in interaction picture with respect to local Hamiltonians, can be
described by the Lindblad form master equation Breuer and Petruccione (2002)

$\displaystyle\dot{\rho}_{AB}$ $\displaystyle=\mathcal{L}_{AB}[\rho_{AB}]$
$\displaystyle=-i[g(ab^{\dagger}+ba^{\dagger})+F\left(e^{-i\Delta_{af}t}a+e^{i\Delta_{af}t}a^{\dagger}\right),\rho_{AB}]$
$\displaystyle+\gamma(N(T)+1)\mathcal{D}_{a}[\rho_{AB}]+\gamma
N(T)\mathcal{D}_{a^{\dagger}}[\rho_{AB}],$ (2)

where $\Delta_{af}=\omega_{f}-\omega$ is the detunning between the laser field
and the charger $A$ and $N(T)$ the average number of photon in mode
$\omega_{a}$ of the environment at temperature $T$ and

$\mathcal{D}_{c}[\rho]=c\rho
c^{\dagger}-\dfrac{1}{2}\\{c^{\dagger}c,\rho_{AB}\\},$ (3)

is the dissipator. In the following it will also be assumed that the charger
and the battery are on resonant with each other
$\omega_{a}=\omega_{b}=\omega$. The results are plotted in Fig. (3) for zero
detuning $\Delta_{af}=0$. As is seen for zero temperature of the environment,
where $N(T)=0$, the maximum energy stored in the battery is almost 0.9, i.e.,
less than one quanta of energy (see Fig. (3a)) and for nonzero temperature of
the environment, with $N(T)=1$, the maximum energy stored in the battery is
around 1.2 (see Fig. (3b)). Therefore, this motivated us to study other
methods of boosting energy transfer into the battery. Those will be explained
in the following sections. We note that for a quantum harmonic oscillator
linearly coupled to a zero-temperature bath, the Markovian dynamics maps
coherent state into coherent state. Since the initial state is chosen to be
vacuum (also a coherent state) the state remains coherent (and hence pure) for
all the times Kossakowski (1972). Hence the ergotropy extraction will be equal
to the internal energy of the battery, i.e., all of the stored energy can be
extracted in the form of work (Fig. (3a)).

Figure 3: (Color online) Energy $E_{B}$ and ergotropy $\mathcal{W}_{B}$ of a
quantum battery $B$, charged by a quantum charger $A$ versus time for a)
$N(T)=0$, $g=0.2$, $F=0.1$, and $\Delta_{af}=0$. A laser field is shining on
the charger and both the charger and the battery are initially in ground
states. Since the state of the battery is pure throughout the process the
ergotropy extraction is equal to the internal energy of the battery. Here it
is assumed that. b) The same as in (a) but for $N(T)=1$.

##  III Tuning the external laser field

Having a thermodynamic system, with an external laser field shining on it, in
which there exists a flow (transfer) of energy the intuitive way to boost this
energy flow (here from the charger to the battery) is to tune the laser field
to be on resonance with the whole system (of the charger-battery). Therefore
introducing the super-modes (global modes) of our charger-battery system

$C_{\pm}=a\pm b,$ (4)

the Hamiltonian in Eq. (II) may be rewritten as

$\displaystyle H$
$\displaystyle=\omega_{+}C_{+}^{\dagger}C_{+}+\omega_{-}C_{-}^{\dagger}C_{-}+\dfrac{F}{\sqrt{2}}\left(e^{-i\omega_{f}t}C_{+}+e^{i\omega_{f}t}C_{+}^{\dagger}\right)$
$\displaystyle+\dfrac{F}{\sqrt{2}}\left(e^{-i\omega_{f}t}C_{-}+e^{i\omega_{f}t}C_{-}^{\dagger}\right),$
(5)

where

$\omega_{\pm}=\omega\pm g.$ (6)

So knowing the coupling strength we can always tune the laser field frequency
to one of these super-modes. In Fig. (4) the energy and ergotropy of the
battery are plotted for cases $\omega_{f}=\omega_{\pm}$ and it is observed
that the ergotropy extraction keeps boundlessly increasing with time. One
issue with Eq. (6) that may not be experimentally friendly is that in order to
find $\omega_{\pm}$ one needs to make a measurement to know the coupling
strength. Introducing the catalyst qubit gives us the chance to avoid this
issue while the same boost in the energy transfer. This will be addressed in
the following section.

Figure 4: (Color online) Energy $E_{B}(t)$ and ergotropy $\mathcal{W}_{B}(t)$
of a quantum battery $B$, charged by a quantum charger $A$ versus time for a)
$\omega_{+}=\omega+g$ (the same plot is obtained for $\omega_{-}=\omega-g$),
$N(T)=0$, $g=0.2$, and $F=0.1$. It is seen that energy transfer in the
resonance case, where $\omega_{f}=\omega_{\pm}$, keeps increasing with time.
Note that since the state of the battery is pure at all times
$E_{B}(t)=\mathcal{W}_{B}(t)$. Both the charger and the battery are initially
in ground states. Since the state of the battery is pure throughout the
process the ergotropy extraction is equal to the internal energy of the
battery. b) The same as in (a) but for $N(T)=1$. In this case energy and
ergotropy are not the same anymore.

##  IV Catalysis in Charging quantum batteries

Figure 5: (Color online) A two-level quantum system interacts with the charger
and the battery, playing the role of a catalyst.

In this section we aim to boost the ergotropy extraction of the battery by
inserting a qubit catalyst in the process of charging between the charger and
the battery (see Fig. (5)). The charger interacts with the qubit and the qubit
interacts with the oscillator while there exists no direct interaction between
the charger and the battery. The Hamiltonian of the whole system becomes

$\displaystyle H$ $\displaystyle=\omega
a^{\dagger}a+\omega_{q}q^{\dagger}q+\omega
b^{\dagger}b+g_{aq}(aq^{\dagger}+a^{\dagger}q)$
$\displaystyle+g_{bq}(bq^{\dagger}+b^{\dagger}q)+F\left(e^{-i\omega_{f}t}a+e^{i\omega_{f}t}a^{\dagger}\right),$
(7)

where $q$ ($q^{\dagger}$) are annihilation (creation) operators of the qubit.
Introducing the super-mode operators

$C_{+}=\sin{\theta}\,a+\cos{\theta}\,b,$ (8)
$C_{-}=\cos{\theta}\,a-\sin{\theta}\,b,$ (9)

where

$\sin{\theta}=\dfrac{g_{aq}}{\sqrt{g_{aq}^{2}+g_{bq}^{2}}},\
\quad\cos{\theta}=\dfrac{g_{bq}}{\sqrt{g_{aq}^{2}+g_{bq}^{2}}},$ (10)

the Hamiltonian in Eq. (II) can be rewritten in the form

$\displaystyle H$
$\displaystyle=\omega_{+}C_{+}^{\dagger}C_{+}+\omega_{-}C_{-}^{\dagger}C_{-}+\omega_{q}q^{\dagger}q+g(C_{+}q^{\dagger}+C_{+}^{\dagger}q)$
$\displaystyle+F\sin{\theta}(e^{-i\omega_{f}t}C_{+}+e^{i\omega_{f}t}C_{+}^{\dagger})$
$\displaystyle+F\cos{\theta}(e^{-i\omega_{f}t}C_{-}+e^{i\omega_{f}t}C_{-}^{\dagger}),$
(11)

where, in this case, for the super-mode frequencies we have
$\omega_{+}=\omega_{-}=\omega$ and $g=\sqrt{g_{aq}^{2}+g_{bq}^{2}}$. Therefore
for zero detuning $\Delta_{af}=\omega-\omega_{f}=0$ the laser field will be on
resonance with the system of charger-battery (see Appendix (VI) for a
generalization to arbitrary number of modes (cells) in the battery). Then a
boost in transfer of energy to the battery will be expected. As illustrated in
Fig. (6), numerical results show that the internal energy and ergotropy
extraction of the battery keep increasing with time with no bound on it. This
is remarkable because here the qubit acts only as a catalyst which facilitates
more transfer of energy from the charger to the battery while it stores almost
no energy during the whole process of charging therefore it can be used for
another charging process. Another interesting point is that the catalyst qubit
does not get entangled with neither the charger and the battery nor the
charger. Of particular note is the operational benefit of the catalyst
approach over the laser-tuning since one does not need to know the coupling
constant in the catalyst approach which makes it more experimentally friendly.
On the hand, it should be noticed that in the dispersive regime,
$\Delta_{aq}=\omega-\omega_{q}\gg g$, where the frequency of the catalyst
qubit is far detuned from the frequency of the charger or the battery energy
transfer decreases and for large detunings it becomes zero.

Figure 6: (Color online) Energies and ergotropy $E_{B}$, $E_{C}$,
$\mathcal{W}_{B}$ of the battery $B$ and the catalyst qubit $Q$, respectively,
versus time for a) $N(T)=0$, $g_{aq}=g_{bq}=0.2$, $F=0.1$, and
$\Delta_{af}=0$. Adding a qubit between the charger and the battery causes the
energy and the ergotropy extraction of the battery to boundlessly keep
increasing with time. As is observed the catalyst qubit does not store any
energy throughout the evolution. Since the state of the battery is pure
throughout the process the ergotropy extraction is equal to the internal
energy of the battery. b) The same as in (a) but for $N(T)=1$.

##  V Summary

We proposed two approaches to boost the energy transfer from a quantum charger
to a quantum battery. The first approach is based on the tuning of the laser
field with the super-mode frequencies of the charger-battery system such that
the laser field will be on resonance with the super-modes of the charger-
battery system. As was shown tuning of the laser-field greatly boosts the
energy transfer in a way that there will be no bound on the work extraction of
the battery. Since the super-modes are dependent on the interaction coupling
constant between the charger and the battery this method requires knowing the
interaction coupling constant. In the second approach a qubit initially in the
ground state, as a catalyst, is added to the system and as a result the laser
field becomes on resonance with the super-modes of the charger-battery system.
Consequently it was seen that the amount of energy transferred to the battery
is greatly boosted again with no bound it, while the qubit interestingly
remains in the ground state throughout almost the whole evolution. As a final
remark, we note that the current proposal can be implemented with the current-
state-of-art superconducting technology Upadhyay et al. (2021). The high level
control of the nanofabrication techniques of the quantum circuits allows
careful deigning and fabrication of qubits and resonators with the coupling
rates reaching strong-coupling regime $g_{aq}\sim\omega_{q}$ Blais et al.
(2021). The possibility to tune the qubit-resonator coupling rate Kafri et al.
(2017) allows to reach the target values needed to see effects described
above. This provides additional control parameter to optimize the energy
transfer between the charger and battery, and ultimately realize the current
platform.

##  VI Appendix

If we have $k$ modes (cells) in the battery and each mode interacts with the
qubit separately the Hamiltonian reads

$\displaystyle H$ $\displaystyle=\sum_{i=0}^{k}\omega
a_{i}^{\dagger}a_{i}+\omega_{q}q^{\dagger}q+\sum_{i=0}^{k}g_{i}(a_{i}q^{\dagger}+a_{i}^{\dagger}q)$
$\displaystyle+F\left(e^{-i\omega_{f}t}a_{0}+e^{i\omega_{f}t}a_{0}^{\dagger}\right).$
(12)

where $a_{0}$ ($a^{\dagger}_{0}$) is the annihilation (creation) operator of
the charger and $a_{i}$ ($a^{\dagger}_{i}$), $i\neq 0$, are the annihilation
(creation) operators of the cells of the battery. Following the Bogoljubov
method Bogoljubov (1958), we define the super-mode operators $C_{i}$ as

$C_{i}=\sum_{l=0,\dots,k}\frac{g_{il}}{n_{l}}a_{l},\ \quad i=0,\dots,k,$ (13)

where $n_{l}=\sqrt{\sum\limits_{j=0,\dots,k}g_{jl}^{2}}$, and $g_{ik}$ are
constants (to be determined), such that the normalized vectors
$\textbf{G}_{i}$ defined by $g_{ik}$ as

$\textbf{G}_{i}=(g_{0i},\ g_{1i},\ \dots,\ g_{ki})/n_{i},\ \quad i=0,\dots,k$
(14)

are mutually orthogonal and

$g_{0l}/n_{l}=g_{l},\ \ l=0,\dots,k.$ (15)

Then the Hamiltonian in Eq. (VI) may be rewritten as

$\displaystyle H$ $\displaystyle=\sum_{i=0}^{k}\omega
C_{i}^{\dagger}C_{i}+\omega_{q}q^{\dagger}q+\left(C_{0}q^{\dagger}+C_{0}^{\dagger}q\right)$
$\displaystyle+F\sum_{i=0}^{k}\chi_{i}\left(e^{-i\omega_{f}t}C_{i}+e^{i\omega_{f}t}C_{i}^{\dagger}\right).$
(16)

in which $\chi$ are functions of coefficients $g_{ik}$. In general, the
transformation in Eq. (13) may be described by a matrix $\hat{G}$ whose $i$-th
column is formed by a vector $\textbf{G}_{i}$ (note that labels of columns
belong to the set $\\{0,1,\dots,k\\}$). We therefore have
$\chi_{i}=[\hat{G}^{-1}]_{0i}$. Thus for $\Delta_{af}=0$ the laser field will
be on resonance with the system of charger-battery.

### Acknowledgements

We acknowledge support from the Foundation for Polish Science through IRAP
project co-financed by EU within the Smart Growth Operational Programme
(contract no.2018/MAB/5). S.B. acknowledges funding by the Natural Sciences
and Engineering Research Council of Canada (NSERC) through its Discovery
Grant, funding and advisory support provided by Alberta Innovates through the
Accelerating Innovations into CarE (AICE) – Concepts Program, and support from
Alberta Innovates and NSERC through Advance Grant.

## References

  * Spohn (1978) H. Spohn, Journal of Mathematical Physics 19, 1227 (1978).
  * Alicki (1979) R. Alicki, Journal of Physics A: Mathematical and General 12, L103 (1979).
  * Parrondo et al. (2015) J. M. Parrondo, J. M. Horowitz, and T. Sagawa, Nature physics 11, 131 (2015).
  * Brandao et al. (2015) F. Brandao, M. Horodecki, N. Ng, J. Oppenheim, and S. Wehner, Proceedings of the National Academy of Sciences 112, 3275 (2015).
  * Ahmadi et al. (2021) B. Ahmadi, S. Salimi, and A. Khorashad, Scientific reports 11, 1 (2021).
  * Ahmadi et al. (2019) B. Ahmadi, S. Salimi, and A. Khorashad, arXiv preprint arXiv:1912.01983 (2019).
  * Alicki and Fannes (2013) R. Alicki and M. Fannes, Physical Review E 87, 042123 (2013).
  * Campaioli et al. (2017) F. Campaioli, F. A. Pollock, F. C. Binder, L. Céleri, J. Goold, S. Vinjanampathy, and K. Modi, Physical Review Letters 118, 150601 (2017).
  * Campaioli et al. (2018) F. Campaioli, F. A. Pollock, and S. Vinjanampathy, _Quantum batteries - review chapter_ (2018), eprint arXiv:1805.05507.
  * Barra (2019) F. Barra, Physical Review Letters 122, 210601 (2019).
  * Farina et al. (2019) D. Farina, G. M. Andolina, A. Mari, M. Polini, and V. Giovannetti, Physical Review B 99, 035421 (2019).
  * Julià-Farré et al. (2020) S. Julià-Farré, T. Salamon, A. Riera, M. N. Bera, and M. Lewenstein, Physical Review Research 2, 023113 (2020).
  * Ferraro et al. (2018) D. Ferraro, M. Campisi, G. M. Andolina, V. Pellegrini, and M. Polini, Physical review letters 120, 117702 (2018).
  * Andolina et al. (2019a) G. M. Andolina, M. Keck, A. Mari, V. Giovannetti, and M. Polini, Physical Review B 99, 205437 (2019a).
  * Andolina et al. (2019b) G. M. Andolina, M. Keck, A. Mari, M. Campisi, V. Giovannetti, and M. Polini, Physical review letters 122, 047702 (2019b).
  * Jonathan and Plenio (1999) D. Jonathan and M. B. Plenio, Physical Review Letters 83, 3566 (1999).
  * Eisert and Wilkens (2000) J. Eisert and M. Wilkens, Physical Review Letters 85, 437 (2000).
  * Loudon and Petruccione (2000) R. Loudon and F. Petruccione, _The Quantum Theory of Light_ (Oxford University Press, 2000).
  * Breuer and Petruccione (2002) H. P. Breuer and F. Petruccione, _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Kossakowski (1972) A. Kossakowski, Reports on Mathematical Physics 3, 247 (1972).
  * Upadhyay et al. (2021) R. Upadhyay, G. Thomas, Y.-C. Chang, D. S. Golubev, A. Guthrie, A. Gubaydullin, J. T. Peltonen, and J. P. Pekola, Physical Review Applied 16, 044045 (2021).
  * Blais et al. (2021) A. Blais, A. L. Grimsmo, S. Girvin, and A. Wallraff, Reviews of Modern Physics 93, 025005 (2021).
  * Kafri et al. (2017) D. Kafri, C. Quintana, Y. Chen, A. Shabani, J. M. Martinis, and H. Neven, Physical Review A 95, 052333 (2017).
  * Bogoljubov (1958) N. N. Bogoljubov, Il Nuovo Cimento 7, 794 (1958).


# A scalable 2-local architecture for quantum annealing of all-to-all Ising
models

Ana Palacios de Luis<EMAIL_ADDRESS>Qilimanjaro Quantum Tech,
Barcelona 08007, Spain  Departament de Física Quàntica i Astrofísica, Facultat
de Física, Universitat de Barcelona, E-08028 Barcelona, Spain  Institut de
Ciències del Cosmos, Universitat de Barcelona, ICCUB, Martí i Franquès 1,
E-08028 Barcelona, Spain.  Artur Garcia-Saez  Qilimanjaro Quantum Tech,
Barcelona 08007, Spain  Barcelona Supercomputing Center, Barcelona 08034,
Spain.  Marta P. Estarellas  Qilimanjaro Quantum Tech, Barcelona 08007, Spain

###### Abstract

Achieving dense connectivities is a challenge for most quantum computing
platforms today, and a particularly crucial one for the case of quantum
annealing applications. In this context, we present a scalable architecture
for quantum annealers defined on a graph of degree $d=3$ and containing
exclusively 2-local interactions to realize an all-to-all connected Ising
model. This amounts to an efficient braiding of logical chains of qubits which
can be derived from a description of the problem in terms of triangles. We
also devise strategies to address the challenges of scalable architectures,
such as the faster shrinking of the gap due to the larger physical Hilbert
space, based on driver Hamiltonians more suited to the symmetries of the
logical solution space. We thus show an alternative route to scale up devices
dedicated to classical optimization tasks within the quantum annealing
paradigm.

††preprint: APS/123-QED

Introduction -  Quantum annealing (QA) [1, 2, 3] is an analog form of quantum
computation which encodes the solution to a certain problem in the ground
state of a final Hamiltonian $H_{f}$. This is done by preparing the ground
state of a trivial Hamiltonian $H_{0}$ and interpolating towards $H_{f}$
following

$H(t)=(1-\lambda(t))H_{0}+\lambda(t)H_{f},$ (1)

where $\lambda(0)=0$ and $\lambda(t_{f})=1$ with $t_{f}$ being the total time
of the anneal.

An interesting Hamiltonian in the context of annealing applications is the
Ising model

$H_{\text{ising}}=H_{f}=\sum_{i}h_{i}\sigma^{z}_{i}+\sum_{j>i}J_{ij}\sigma^{z}_{i}\sigma^{z}_{j},$
(2)

which is relatively simple to implement in a number of quantum computing
platforms, such as superconducting quantum circuits [4, 5] or neutral atoms
[6, 7]. The cost function of many relevant classical problems can be cast in
the shape of (2), often requiring all-to-all interactions. However,
establishing a dense connectivity among all qubits remains a challenging task,
since the hard-wiring of all the necessary links is not scalable due to
crosstalk and packing issues within the chip. Based on [8], it was proposed in
[9, 10] to make up for the required connectivity by establishing chains of
qubits that represent the same logical variable. This is then implemented in a
fixed family of graphs within which one searches for an optimal embedding, an
approach known as minor embedding. Another alternative [11] is switching to
the description of the problem in terms of parities between the original
variables. This allows to encode the original cost function into single-qubit
terms at the expense of additional 4- and 3-body interactions (or 2-local
interactions mediated by qutrits), which are harder to achieve experimentally
(however, some proposals do exist [12]). Recently, a new method based on
perturbative gadgets has been proposed for realizing all-to-all connectivity
[13], but which comes at the cost of interaction strengths of
$\mathcal{O}(N^{6})$ for an $N$-qubit problem.

In this letter, we provide a novel architecture that enables the encoding of
arbitrary Ising models, while requiring only 2-local interactions, where every
qubit is linked to other 3 qubits at most. We develop this architecture step
by step from a decomposition of the original problem in terms of triads of
qubits, and in doing so we naturally arrive at a family of alternative
formulations of the original problem of independent interest, which were
previously pointed out in [14, 15].

Triangle decomposition of the Ising problem -  In order to kick-start the
description of our encoding organically, we first consider the less general
class of $\mathbb{Z}_{2}$-symmetric classical Ising models for a system of
size $N$, described as:

$\displaystyle
H_{\mathbb{Z}_{2}}=\sum_{i,j}J_{ij}\sigma^{z}_{i}\sigma^{z}_{j}$ (3)

The starting point of this construction is the fact that the energy of an all-
to-all connected Hamiltonian of the shape (3) can be described as the sum of
the energies of single triangles with edges $j^{\triangle}_{ij}$, as long as
$J_{ij}=\sum_{\triangle\supset\\{i,j\\}}j^{\triangle}_{ij}$.

$\displaystyle H_{\mathbb{Z}_{2}}=\sum_{\triangle}H_{\triangle}$ (4)
$\displaystyle
H_{\triangle}=j^{\triangle}_{ij}\sigma^{z}_{i}\sigma^{z}_{j}+j^{\triangle}_{ik}\sigma^{z}_{i}\sigma^{z}_{k}+j^{\triangle}_{jk}\sigma^{z}_{j}\sigma^{z}_{k}\text{
for }\\{i,j,k\\}\subset\triangle$ (5)

$q_{0}q_{1}$ | $00$ | $01$ | $10$ | $11$
---|---|---|---|---
$l_{0}l_{1}l_{2}$ |  $000$, $111$ |  $110$, $001$ |  $011$, $100$ |  $101$, $010$

$q_{0}q_{1}$ | $00$ | $01$ | $10$ | $11$
---|---|---|---|---
$q_{2}$ | $0$ | $1$ | $0$ | $1$ | $0$ | $1$ | $0$ | $1$
$l_{0}l_{1}l_{2}$ | $000$ | $111$ | $110$ | $001$ | $011$ | $100$ | $101$ | $010$

Table 1: Encoding of logical triangle configurations ($l_{i}$ variables) in
physical spins ($q_{i}$) for the $\mathbb{Z}_{2}$-symmetric case (a) and non-
symmetric case (b) including local fields, which includes an additional sign
qubit $q_{2}$.

To describe the encoding of triangular cells we employ 2 qubits, since this is
enough to describe the energetically distinct states of the model due to
inversion symmetry. This is illustrated in Table 1a. a)b)

Without loss of generality, we adopt the convention of always labeling the
all-parallel configuration as “00”. In the following description, we refer to
the space in which we are representing triangles in terms of pairs of spins as
the physical representation, while the original model is referred to as the
logical representation.

We can determine the Hamiltonian $H_{\triangle}^{\text{physical}}$ acting on a
physical qubit pair that reproduces the correct energies by solving a simple
$4\times 4$ linear system. Following Table 1a. we find:

$\displaystyle
H_{\triangle}^{\text{physical}}=h_{0}\sigma^{z}_{q_{0}}+h_{1}\sigma^{z}_{q_{1}}+J\sigma^{z}_{q_{0}}\sigma^{z}_{q_{1}}$
(6) $\displaystyle\text{where}\qquad h_{0}=j_{01},\qquad h_{1}=j_{02},\qquad
J=j_{12}.$ (7)

Consistency conditions among triangles -  In order to ensure that the separate
triangular cells refer to a logical state in the original model, we enforce
some constraints. These can be deduced from the realization that shared edges
among triangles, also referred to as interacting edges from here on, must be
consistent with each other: if the configuration of one cell is such that the
shared edge holds parallel spins, this must also be the case for other
triangular cells that contain said edge. It is easy to check that, for the
encoding illustrated in Table. 1a., this can always be achieved by the means
of strong 2-local ferromagnetic interactions, as long as every cell has at
least one edge that is not shared with any other cell. This non-interacting
edge is labeled with the “11” configuration in order to avoid the need for
many-body terms. With this, the full physical Hamiltonian embedding the
logical problem is:

$\displaystyle
H_{p}=\sum_{\triangle}H_{\triangle}^{\text{physical}}+H_{\text{constraints}}$
(8) $\displaystyle
H_{\text{constraints}}=\sum_{p\in\\{\triangle\\}}\sum_{\begin{subarray}{c}q\in\\{\triangle\\}\\\
q>p\end{subarray}}J_{P}(I-\sigma^{z}_{p_{i}}\sigma^{z}_{q_{j}})$ (9)

where $i,j$ can be 0 or 1 depending on the internal labeling of the triangular
cell configurations. For sufficiently strong ferromagnetic penalties $J_{P}$,
the Hamiltonian (8) exactly encodes the spectrum of the original
$H_{\mathbb{Z}_{2}}$ in its low-energy eigenspace. We note that the generation
of strong enough ferromagnetic interactions as system size grows is an issue
that presents itself in all penalty-enforced embedding. Some studies on the
scalability of penalties have been carried out for the minor embedding and LHZ
cases, showing linear [16, 17] and up to quadratic [17] penalty growth in the
range of problems and system sizes examined; similar results are expected for
the present scheme due to its similarities with them. We tackle this
limitation further on in this paper by proposing driver Hamiltonians that ease
these penalty strength requirements.

Scaling up the construction -  For an arbitrary graph of size $N$ we can
always find some family of decompositions in terms of triangles such that
every edge is accounted for at least once, and such that all triangles hold a
non-interacting edge. This is achieved by choosing all the triangles such that
they have a selected node $k^{*}$ in common. Starting from a single triangular
cell, the graph can be increased from size $\eta$ to $\eta+1$ according to
this criterion as follows: for the addition of each new node $i=\eta+1$, we
include the triangles that complete the missing links to restore full
connectivity by selecting the triads $\\{i,k,k^{*}\\}$, where
$k=1,...,k^{*}-1,k^{*}+1,...,\eta$ is one of the remaining older nodes. With
this scheme, the addition of an $\eta$-th qubit adds $\eta-2$ new triangles,
and thus $2(\eta-2)$ physical qubits, such that the total number of qubits $n$
required by the physical implementation is:

$n=2\sum_{\eta=1}^{N-2}\eta=(N-1)(N-2)$ (10)

To complete the algorithmic formulation of the construction, we consider the
inclusion of the full coupling strength between two qubits in the first
triangular cell where it appears, and set it to zero in the following
triangular cells that contain said edge. We note that this choice is valid
because it satisfies the condition
$J_{ij}=\sum_{\triangle}j^{\triangle}_{ij}$. With this, the procedure to build
the full hardware graph can be rephrased as follows. First, we choose the
selected node $k^{*}$, i.e., the logical variable to be involved in all the
triangles. Then, we arrange the pairs of physical qubits representing
triangular cells as follows: the first (top-most) row contains the cell
$(0,1,k^{*})$, with $k^{*}$ being the selected node, the second row contains
$(0,2,k^{*})$ and $(1,2,k^{*})$, etc., skipping over $k^{*}$. The resulting
graph is illustrated in Fig. 1 for the case of $N=5$ and $k^{*}=2$. Following
the consistency conditions previously described, we end up with strong
ferromagnetic links connecting neighboring cells that share an edge in such a
way that every physical qubit is connected to at most other 3, taking into
account the link that connects the pair of physical qubits within each cell.

Incorporation of local fields -  We can simulate a general Ising model like
the one in Eq. (2) by slightly modifying the current construction: due to the
symmetry breaking, now an additional physical qubit is required for the
description of each triangle. We refer to the latter as a sign qubit, and it
indicates whether the parallel spins are pointing upwards or downwards
(represented by $q_{2}$ in Table 1b.). The single-cell Hamiltonian giving the
right energy contributions can be found by solving the $8\times 8$ linear
system, in accordance with what was done in the $\mathbb{Z}_{2}$-symmetric
case.

FIG. 1: Full hardware graph representing a sparse realization of Eq. (2) for
$N=5$, which also corresponds to a sparse realization of Eq. (3) for $N=5$ in
the absence of the sign qubits. Blue dotted circles group the qubits
representing a single cell, corresponding to the triad of qubits in blue
closest to it. The purple dashed links correspond to strong ferromagnetic
connections among sign qubits. Orange and (dashed) green color-coding
corresponds to the $\\{J_{ij}\\}$ and $\\{h_{i}\\}$ of the logical problem. In
this illustration, $k^{*}=2$.

The consistency conditions regarding sign qubits imply that they are
ferromagnetically coupled to the sign qubits of their neighboring cells.
However, following again the rule of including the full local field on the
first appearance, the only sign qubits that maintain some connectivity with
the rest of their cell are those corresponding to the left-most diagonal.
Thus, the remaining sign qubits are expendable and can be spared by collapsing
ferromagnetic links. The resulting scheme is illustrated for $N=5$ in Fig. 1.

We highlight that if the logical chains are allowed to be of length $N-1$
instead of $N-2$, the incorporation of local fields can be done within the
same structure as the one described in the absence of local fields.

Driver Hamiltonians for scalability -  Two issues arise when embedding a
problem within a larger Hilbert space by means of some constraints: the
preservation of the logical ground state as the ground state of $H_{p}$ (which
amounts to enforcing the constraints strongly enough) and the scaling of the
minimal gap, since it now decreases with the physical system size $n$ rather
than the original logical system size $N$. The standard driver Hamiltonian for
annealing $H_{0}=H_{std}=-\sum_{i}\sigma^{x}_{i}$ encodes no information about
the feasible solution subspace, and is thus expected to render a search
through a space much larger than necessary. This makes the computation less
efficient and translates into an exponentially decreasing minimal gap in
comparison to the fully connected model. This was also pointed out in the SQA
study benchmarking the performance of minor embedding vs. LHZ in [18].

However, we do have information about the relevant subspace to explore: it is
the space spanned by $\\{|0\rangle^{\otimes l},|1\rangle^{\otimes l}\\}$
within each chain, which we will refer to as the logical subspace
$\mathcal{L}$. The collection of optimal driver Hamiltonians for the annealing
of $H_{p}$ 111 In the sense of restricting the evolution to the feasible
space.  should thus commute with the operator $P$ describing the restriction
to this subspace [20, 21], i.e., $[P,H_{0}]=[P,H_{p}]=0$, $[H_{0},H_{p}]\neq
0$ and also induce hopping between all regions of the solution space. Notice
that for our architecture $P=\sum_{c_{k}\in\mathcal{C}}\sum_{i\in
c_{k}}\sigma^{z}_{i}\sigma^{z}_{i+1}$. These conditions automatically solve
both of the previously described issues by removing the need for constraints
and making the annealing of the embedded problem equivalent to that of the
fully connected model. However, a Hamiltonian of these characteristics for the
class of $H_{p}$ at hand will have highly nonlocal terms 222 Exemplary
Hamiltonians of this type can be extracted from a circuit that builds a GHZ
state. One such Hamiltonian is:  $\displaystyle
H_{0}=H_{\text{GHZ}}=-\sum_{c_{k}\in\mathcal{C}}\left(\prod_{i_{r}\in
c_{k}}\sigma^{x}_{i_{r}}+\sigma^{z}_{i_{0}}\sum_{r=1}^{l-1}\sigma^{z}_{i_{r}}\right)$
(11) where $c_{k}$ corresponds to the $k$-th chain of the chain set
$\mathcal{C}$, all of which are of length $l$. The ground and first excited
states of Eq. (16) correspond to the even and odd GHZ states. , rendering its
physical implementation impractical.

In order to find more feasible alternatives, we relax the previous conditions
such that the commutation of the driver with the constraint is no longer
satisfied, but such that the low-energy subspace of the driver has a higher
overlap with the logical subspace $\mathcal{L}$.

The easiest driver to implement that would fulfill the previous requirement is
the ferromagnetic transverse-field Ising model (TFIM):

$\displaystyle
H_{0}=H_{\text{TFIM}}=\sum_{c_{k}\in\mathcal{C}}h_{\text{TFIM}}^{c_{k}}$ (12)
$\displaystyle
h_{\text{TFIM}}^{c_{k}}=-\sum_{i}\sigma^{x}_{i}-J_{\text{ZZ}}\sum_{i\in
c_{k}}\sigma^{z}_{i}\sigma^{z}_{i+1}$ (13)

In the regime where the local field can be considered as a perturbation to the
interactions, $J_{\text{ZZ}}\gg 1$, it is well known that the two lowest
excited states approach the symmetric and antisymmetric GHZ states (see [23],
for example). However, in this regime the gap between first and ground states
also vanishes exponentially. Thus, in practice we can only approach the target
subspace $\mathcal{L}$ from the gapped region $J_{\text{ZZ}}<1$, where the gap
tends to $2(1-J_{\text{ZZ}})$ asymptotically, and get as close as possible to
the transition to the ferromagnetic phase at $J_{\text{ZZ}}=1$, which is
gapless in the thermodynamic limit [24].

A better scaling alternative comes with the inclusion of additional
interactions in the driver Hamiltonian:

$\displaystyle
H_{0}=H_{\text{XYZ}}=\sum_{c_{k}\in\mathcal{C}}h_{\text{XYZ}}^{c_{k}}$ (14)
$\displaystyle
h_{\text{XYZ}}^{c_{k}}=-\sum_{i}\sigma^{x}_{i}-J_{\text{ZZ}}\sum_{i\in
c_{k}}(\sigma^{z}_{i}\sigma^{z}_{i+1}+\alpha\sigma^{x}_{i}\sigma^{x}_{i+1})$
(15)

The previous driver is nothing but the XYZ model, which is integrable in the
absence of local fields [25]. In this latter case, it has been shown to
present a ferromagnetic phase in the $z$-direction for $\alpha<\frac{1}{2}$
[26]. We numerically study the model in (15) with DMRG in order to evaluate
the scaling of the gap and of the projection onto $\mathcal{L}$ of the ground
state $\langle{P_{\mathcal{L}}}\rangle_{GS}$, with
$P_{\mathcal{L}}=|0\rangle^{\otimes l}\langle 0|^{\otimes
l}+|1\rangle^{\otimes l}\langle 1|^{\otimes l}$ and $l$ the length of the
chain. The intuition behind this choice of Hamiltonian is that the presence of
the XX term will enlarge the gap and thus offer more favorable routes of
approaching the gapless ferromagnatic phase. The decay of the ground state
overlap with $\mathcal{L}$ can be fitted to an exponential of the form
$\langle GS|P_{\mathcal{L}}|GS\rangle=2^{-\beta l+c}$. Notice that for
$H_{std}$, $\beta=1$ and $c=1$. Thus, a smaller $\beta$ exponent indicates a
bias towards the logical subspace; the smaller $\beta$, the greater the bias.
Fig. 2 shows how appropriate combinations ($\alpha,J_{\text{ZZ}}$) can achieve
a reduction of $\beta$ of around an order of magnitude. The comparison in Fig.
2 also manifests the greater flexibility of the XYZ Hamiltonian with respect
to the TFIM, evidenced in the fact that tuning the relative strength of the XX
couplings allows for achieving a greater bias towards $\mathcal{L}$ for equal
gap sizes. The presence of these additional couplings also allows to expand
the range of $J_{\text{ZZ}}$ to be implemented in comparison to that of the
TFIM, such that better biases can be achieved for a still acceptable gap size.
In this manner, we can achieve a reduction of about an order of magnitude in
$\beta$.

FIG. 2: Decay coefficient of the projection onto $\mathcal{L}$, $\beta$, vs.
the gap $\Delta$ for different values of $J_{\text{ZZ}}$. Results are relative
to a single chain of $l=300$, where it was numerically checked that the gap
had stabilized to its value in the thermodynamic limit.
$J_{\text{ZZ}}\in[2,4]$ for $h_{\text{XYZ}}$ and $J_{\text{ZZ}}\in[0.8,0.98]$
for $h_{\text{TFIM}}$, always in equal intervals, with smaller gaps
corresponding to higher $J_{\text{ZZ}}$.

Discussion -  We have obtained a simple hardware graph of maximum degree $d=3$
and 2-local interactions that is sufficient to encode any arbitrary Ising
problem. By building the architecture from a decomposition of the original
problem, a dual interpretation of the resulting graph is achieved: one in
terms of consistent triangular cells and one in terms of chains of logical
variables. The full architecture, incorporating local fields, requires
$n=N(N-2)$ physical qubits and $c=N(3N-5)/2$ couplers, out of which
$c_{\text{FM}}=N(N-2)$ will anneal to a strong ferromagnetic coupling. The
sparsity of the graph also implies less crosstalk to keep under control, an
important scalability feature. Achieving strong couplings without compromising
coherence is a challenging task, but one that needs to be tackled nonetheless
for any embedding strategy that enables scalability. In the particular case of
superconducting technology, already standard couplers naturally provide a much
larger ferromagnetic range than the antiferromagnetic one [4]. We have
contributed to this scalability challenge from the algorithmic side by
proposing alternative driver Hamiltonians that enhance the baseline success
probability of the anneal by directing the search towards the logical
subspace, mitigating the need for strong penalties. Some of these Hamiltonians
can be readily implemented with current devices.

The triangle architecture discussed here is minimal: a graph of degree $d=3$
is the smallest to hold NP-complete problems [27]. As a consequence, a faulty
element will truncate the maximum embedding size sharply depending on where
the fault is located. Thus, this architecture has no tolerance to low
fabrication yield, which limits its large-scale implementation as of today. In
the particular case of superconducting quantum technologies (specifically for
flux qubits) robust fabrication processes are still in their development
phase, but yield should stop being an issue once the technology is mature.

All in all, the proposed scheme has several advantages: the embedding of the
problem does not require any additional computation, the sparsity of the graph
aids noise and crosstalk minimization and the interactions required are within
the most employed model in current quantum annealing devices. In addition, the
repetition code present in this architecture allows to spatially spread the
information, making the computation more robust to some external error
sources. Finally, the structure of the logical space allows to devise
relatively simple driver Hamiltonians that mitigate the strength of the
constraint terms.

Acknowledgments -  We thank Ramiro Sagastizabal, Matthias Werner, Arnau Riera
and Josep Bosch for valuable discussions. This work was supported by European
Commission EIC-Transition project RoCCQeT (GA 101112839), the Agencia de
Gestió d’Ajuts Universitaris i de Recerca through the DI grant (No. DI74) and
the Spanish Ministry of Science and Innovation through the DI grant (No.
DIN2020-011168).

## References

  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse Ising model, Physical Review E 58, 5355 (1998).
  * Farhi  _et al._ [2000] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, Quantum Computation by Adiabatic Evolution,  (2000), arXiv:quant-ph/0001106.
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Adiabatic quantum computation, Rev. Mod. Phys. 90, 64 (2018).
  * Weber  _et al._ [2017] S. J. Weber, G. O. Samach, D. Hover, S. Gustavsson, D. K. Kim, A. Melville, D. Rosenberg, A. P. Sears, F. Yan, J. L. Yoder, W. D. Oliver, and A. J. Kerman, Coherent coupled qubits for quantum annealing, Physical Review Applied 8, 014004 (2017).
  * Krantz  _et al._ [2019] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, A Quantum Engineer’s Guide to Superconducting Qubits, Applied Physics Reviews 6, 021318 (2019), arXiv: 1904.06560.
  * Pohl  _et al._ [2010] T. Pohl, E. Demler, and M. D. Lukin, Dynamical Crystallization in the Dipole Blockade of Ultracold Atoms, Physical Review Letters 104, 043002 (2010).
  * de Léséleuc  _et al._ [2018] S. de Léséleuc, S. Weber, V. Lienhard, D. Barredo, H. P. Büchler, T. Lahaye, and A. Browaeys, Accurate Mapping of Multilevel Rydberg Atoms on Interacting Spin-$1/2$ Particles for the Quantum Simulation of Ising Models, Physical Review Letters 120, 113602 (2018).
  * Kaminsky  _et al._ [2004] W. M. Kaminsky, S. Lloyd, and T. P. Orlando, Scalable Superconducting Architecture for Adiabatic Quantum Computation,  (2004), arXiv:0403090.
  * Choi [2008] V. Choi, Minor-embedding in adiabatic quantum computation: I. The parameter setting problem, Quantum Information Processing 7, 193 (2008).
  * Choi [2011] V. Choi, Minor-embedding in adiabatic quantum computation: II. Minor-universal graph design, Quantum Information Processing 10, 343 (2011).
  * Lechner  _et al._ [2015] W. Lechner, P. Hauke, and P. Zoller, A quantum annealing architecture with all-to-all connectivity from local interactions, Science Advances 1, e1500838 (2015).
  * Puri  _et al._ [2017] S. Puri, C. K. Andersen, A. L. Grimsmo, and A. Blais, Quantum annealing with all-to-all connected nonlinear oscillators, Nature Communications 8, 15785 (2017).
  * Mozgunov [2023] E. Mozgunov, Precision of quantum simulation of all-to-all coupling in a local architecture,  (2023), arXiv:2302.02458.
  * Fujii  _et al._ [2022] T. Fujii, K. Komuro, Y. Okudaira, R. Narita, and M. Sawada, Energy landscape transformation of Ising problem with invariant eigenvalues for quantum annealing,  (2022), arXiv:2202.05927.
  * Fujii  _et al._ [2023] T. Fujii, K. Komuro, Y. Okudaira, and M. Sawada, Eigenvalue-Invariant Transformation of Ising Problem for Anti-Crossing Mitigation in Quantum Annealing, Journal of the Physical Society of Japan 92, 044001 (2023).
  * 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, Physical Review X 5, 031040 (2015).
  * Lanthaler and Lechner [2021] M. Lanthaler and W. Lechner, Minimal constraints in the parity formulation of optimization problems, New Journal of Physics 23, 083039 (2021).
  * Albash  _et al._ [2016] T. Albash, W. Vinci, and D. A. Lidar, Simulated-quantum-annealing comparison between all-to-all connectivity schemes, Physical Review A 94, 022327 (2016).
  * Note [1] In the sense of restricting the evolution to the feasible space.
  * Hen and Spedalieri [2016] I. Hen and F. M. Spedalieri, Quantum Annealing for Constrained Optimization, Physical Review Applied 5, 034007 (2016).
  * Hen and Sarandy [2016] I. Hen and M. S. Sarandy, Driver Hamiltonians for constrained optimization in quantum annealing, Physical Review A 93, 062312 (2016).
  * Note [2] Exemplary Hamiltonians of this type can be extracted from a circuit that builds a GHZ state. One such Hamiltonian is:  $\displaystyle H_{0}=H_{\text{GHZ}}=-\sum@\slimits@_{c_{k}\in\mathcal{C}}\left(\prod@\slimits@_{i_{r}\in c_{k}}\sigma^{x}_{i_{r}}+\sigma^{z}_{i_{0}}\sum@\slimits@_{r=1}^{l-1}\sigma^{z}_{i_{r}}\right)$ (16) where $c_{k}$ corresponds to the $k$-th chain of the chain set $\mathcal{C}$, all of which are of length $l$. The ground and first excited states of Eq. (16) correspond to the even and odd GHZ states.
  * Mbeng  _et al._ [2020] G. B. Mbeng, A. Russomanno, and G. E. Santoro,  _The quantum Ising chain for beginners_, Tech. Rep. arXiv:2009.09208 (arXiv, 2020) arXiv:2009.09208 [cond-mat, physics:quant-ph] type: article.
  * Pfeuty [1970] P. Pfeuty, The one-dimensional Ising model with a transverse field, Annals of Physics 57, 79 (1970).
  * Baxter [1990] R. J. Baxter,  _Exactly solved models in statistical mechanics_ , 2nd ed. (Acad. Press, London, 1990).
  * Shi  _et al._ [2020] Q.-Q. Shi, S.-H. Li, and H.-Q. Zhou, Duality and ground-state phase diagram for the quantum XYZ model with arbitrary spin $s$ in one spatial dimension, Journal of Physics A: Mathematical and Theoretical 53, 155301 (2020).
  * Barahona [1982] F. Barahona, On the computational complexity of Ising spin glass models, Journal of Physics A: Mathematical and General 15, 3241 (1982).