# Common envelope evolution and triple dynamics as potential pathways to form
the inner white dwarf + brown dwarf binary of the triple star system Gaia
0007$-$1605

Felipe Lagos1, Monica Zorotovic2, Matthias R. Schreiber3,4 B. T. Gänsicke,1,5
1 Department of physics, University of Warwick, Gibbet Hill, Coventry CV4 7AL,
UK
2 Instituto de Física y Astronomía, Universidad de Valparaíso, Av. Gran
Bretaña 1111, Valparaíso, Chile
3 Departamento de Física, Universidad Técnica Federico Santa María, Avenida
España 1680, Valparaíso, Chile
4 Millennium Nucleus for Planet Formation, NPF, Valparaíso, 2340000, Chile
5 Centre for Exoplanets and Habitability, University of Warwick, Coventry CV4
7AL, UK
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The recently discovered system Gaia 0007$-$1605 consisting of a white dwarf
with a close brown dwarf companion and a distant white dwarf tertiary very
much resembles the triple system containing the first transiting planet
candidate around a white dwarf ever discovered: WD 1856+534. We have
previously argued that the inner binary in WD 1856+534 most likely formed
through common envelope evolution but triple star dynamics represent an
alternative scenario. Here we analyze different formation scenarios for Gaia
0007$-$1605\. We reconstructed the potential common envelope evolution of the
system and find that assuming standard parameters for the energy budget
provides a reasonable solution. In agreement with other close white dwarf $+$
brown dwarf binaries, and in contrast to WD 1856+534, no energy sources other
than orbital energy during common envelope evolution are required to
understand the current configuration of the system. In addition, using
analytical prescriptions for triple dynamics, we show that Von
Zeipel–Lidov–Kozai oscillations might have trigger tidal migration due to high
eccentricity incursions ($e\gtrsim 0.997$). We conclude that the inner binary
in Gaia 0007$-$1605, as its sibling WD 1856+534, formed either through common
envelope evolution, triple dynamics or a combination of both mechanisms.

###### keywords:

binaries: close – white dwarfs – brown dwarfs

††pubyear: 2022††pagerange: Common envelope evolution and triple dynamics as
potential pathways to form the inner white dwarf + brown dwarf binary of the
triple star system Gaia 0007$-$1605–Common envelope evolution and triple
dynamics as potential pathways to form the inner white dwarf + brown dwarf
binary of the triple star system Gaia 0007$-$1605

##  1 Introduction

White dwarfs that are members of close binaries111through the paper we use the
term “close binary” to refer to binaries with orbital periods up to hundred of
days. are important in a wide range of astrophysical contexts including
studies of supernovae Ia or the detection of gravitational waves. The
classical formation scenario for close binary stars containing a white dwarf
is common envelope evolution (Paczynski, 1976).

Common envelope evolution is an inherently complicated process and it has so
far turned out to be impossible to cover the large range of spatial and
temporal scales in hydrodynamic simulations (Ivanova et al., 2013). To compare
observations and model predictions, proper simulations of the process are
therefore often replaced by a parameterized energy equation relating the
binding energy of the envelope and the change in orbital energy of the binary.

Observed samples of close white dwarf binaries with main sequence star
companions can be used to constrain the energy budget during common envelope
evolution. In the vast majority of cases, observed populations can be
understood using a small common envelope efficiency and without assuming
contributions from additional energy sources such as recombination energy.
This finding holds for post common envelope binaries with M-dwarf (Zorotovic
et al., 2010), substellar (Zorotovic & Schreiber, 2022), as well as sun-like
(Hernandez et al., 2021, 2022a, 2022b) companions.

However, the general success of the common envelope scenario in explaining the
observed populations of close binaries does not exclude the existence of
alternative formation scenarios of close white dwarf binaries. Recently, Lagos
et al. (2022) showed that the period distribution of white dwarfs with close
G-type companions can be explained if systems with periods of months to years
form through stable but non-conservative mass transfer (see also Masuda et
al., 2019). The situation is similar for the observed population of double
white dwarfs, i.e. their characteristics can be best explained by considering
the possibility of common envelope evolution and stable non-conservative mass
transfer (e.g. Nelemans et al., 2000; Webbink, 2008; Woods et al., 2012).

Another alternative for forming close white dwarf binaries is triple dynamics.
The inner binary of a hierarchical triple system may exchange angular momentum
with the orbit of the distant tertiary through the so-called Von Zeipel-Lidov-
Kozai (ZLK) mechanism (Von Zeipel, 1910; Lidov, 1962; Kozai, 1962), which may
generate large eccentricities and subsequent tidal decay in the inner binary.
ZLK oscillations have been used to explain a large variety of phenomena,
including the formation of hot Jupiters (e.g. Wu & Murray, 2003; Naoz et al.,
2011; Petrovich, 2015) or blue stragglers (e.g. Perets & Fabrycky, 2009) and
may even result in the merger of double white dwarf binaries (e.g. Thompson,
2011).

So far, however, the formation of close white dwarf binaries through ZLK
oscillations remains hypothetical. To the best of our knowledge, a triple
system where an inner binary containing a white dwarf has formed most likely
through triple dynamics and tidal migration has not yet been identified. The
perhaps most promising candidate known is the transiting gas giant planet
around WD 1856+534 which is part of a hierarchical triple system (Vanderburg
et al., 2020). While ZLK oscillations could in principle be responsible for
the currently observed configuration (e.g. Muñoz & Petrovich, 2020; O’Connor
et al., 2021; Stephan et al., 2021), it has also been shown that common
envelope evolution can explain its tight orbit (Lagos et al., 2021; Chamandy
et al., 2021; Merlov et al., 2021). Measuring the mass of the transiting
planet could constrain the evolutionary history of the system. Whilst common-
envelope evolution is disfavoured for masses below $5$ Mjup ZLK oscillations
also cover lower masses being the most likely scenario below $3$ Mjup.
Unfortunately, according to the lower limits derived from transmission
spectroscopy that appear in the literature ($0.84$ and $2.4$ Mjup at the
$2\sigma$ level, Alonso et al. 2021; Xu et al. 2021) both scenarios remain
plausible.

Recently, a system very similar to WD 1856+534 has been found and
characterized by Rebassa-Mansergas et al. (2022). They showed that the
infrared excess of the white dwarf Gaia 0007$-$1605 is caused by a brown dwarf
companion and using spectral fits combined with Gaia photometry they
convincingly constrained most of the parameters of both components, with
exception of the brown dwarf mass and the total age of the system (see Table
1). El-Badry et al. (2021) and Rebassa-Mansergas et al. (2021) furthermore
revealed that the close white dwarf plus brown dwarf (WD+BD) binary is in fact
the inner binary of a hierarchical triple system by discovering a common
proper motion white dwarf companion. The distant white dwarf was characterized
based on photometry only and the estimated parameters – although in agreement
with previous rough estimates (Gentile Fusillo et al., 2021) – are therefore
less certain. The total age of the system derived from the parameters of the
tertiary is $\sim 10$ Gyr which is in agreement with the likely membership of
the system to the Galactic disc.

We here present an investigation of the potential evolutionary history of Gaia
0007$-$1605 using standard prescriptions for common envelope evolution and
analytical approximations describing possible ZLK oscillations. We find that
common envelope evolution represents a natural explanation for the current
configuration and that the formation of Gaia 0007$-$1605 can be reproduced
with a small common envelope efficiency, such as virtually all post common
envelope binaries. Tidal migration induced by ZLK oscillations also represents
a plausible scenario as long as the brown dwarf survives the post-main-
sequence evolution of its host star and its eccentricity reaches values
greater than $\approx 0.997$.

Table 1: Stellar parameters and ages for Gaia 0007$-$1605 reported by Rebassa-
Mansergas et al. (2022). Values with “$*$” require further confirmation. In
particular, the total age of the systems is based on the sum of the cooling
age of the outer white dwarf and its main sequence progenitor lifetime, which
in turn depends on the initial-to-final mass relation. However, the latter is
not well defined in the mass range of the outer white dwarf. The mass of the
brown dwarf is based on its L$3\pm 1$ spectral type.

Parameter | Inner WD | Brown dwarf | Outer WD
---|---|---|---
Mass [M⊙] | $0.54\pm 0.01$ | $0.07^{*}$ | $0.56\pm 0.05$
Orbital period [d] | $1.0446\pm 0.0015$ | $1.0446\pm 0.0015$ | -
Cooling time [Gyr] | $0.360\pm 0.002$ | - | $8.2\pm 0.2$
Total age [Gyr] | $\simeq 10^{*}$ | $\simeq 10^{*}$ | $\simeq 10^{*}$
Projected separation [au] | - | - | 1673.11

##  2 Reconstructing the evolution of the inner binary with common envelope
evolution

We reconstructed the evolutionary history of the inner WD+BD binary in Gaia
0007$-$1605 assuming that it evolved to its current short orbital period
during a common-envelope phase and without influence of the distant companion
on the orbital period decrease. We used the same method recently described in
Zorotovic & Schreiber (2022) for a sample of well characterized close WD+BD
binaries. In summary, we first calculated the period the system had
immediately after ejecting the envelope based on the white dwarf cooling age
and angular momentum loss in the post-common-envelope phase through
gravitational radiation only (Schreiber & Gänsicke, 2003). We then searched
for possible progenitors of the white dwarf within a grid calculated with the
single-star-evolution (SSE) code from Hurley et al. (2000) for solar
metallicity, and reconstructed the common-envelope phase using Roche-geometry
and the energy formalism developed by Webbink (1984), assuming that no energy
sources other than orbital and thermal energy contributed to unbinding the
envelope (for more details see Zorotovic & Schreiber 2022).

In Fig. 1 we present the results for the estimated total age of the system
(top panel), initial mass of the progenitor of the white dwarf
($\mathrm{M}_{\mathrm{prog,wd}}$, middle panel), and orbital period at the
onset of the common-envelope phase ($\mathrm{P}_{\mathrm{CE,in}}$, bottom), as
a function of the common-envelope efficiency $\alpha_{\mathrm{CE}}$ (i.e. the
fraction of the change in orbital energy that is used to unbind the envelope).
The darker results correspond to solutions assuming a white dwarf mass within
the error given by Rebassa-Mansergas et al. (2022, i.e.
$0.01\,\mbox{M${}_{\odot}$}$), while the results marked using light gray allow
the white dwarf mass to vary within a range of $\pm
0.05\,\mbox{M${}_{\odot}$}$. In all cases the white dwarf progenitor filled
its Roche Lobe on the asymptotic giant branch (AGB), meaning the white dwarf
is composed of carbon and oxygen. Assuming that the white dwarf mass is
accurate within the small error given by Rebassa-Mansergas et al. (2022) we
derive a maximum age of $\sim 6$ Gyr for the system, which is not consistent
with the large age the authors derived based on the distant white dwarf
companion. In order to obtain a larger total age of $\sim 10$ Gyr the white
dwarf mass needs to be slightly smaller than estimated by Rebassa-Mansergas et
al. (2021, $\sim 0.51$ M⊙), implying it descends from a low-mass progenitor
($\raisebox{-1.72218pt}{$\stackrel{{\scriptstyle<}}{{\scriptstyle\sim}}$}1.1$
M⊙) with a larger main-sequence lifetime. Also, a low-mass progenitor would
have less mass in the envelope at the onset of the common envelope phase,
which allows the common-envelope efficiency to be smaller, i.e.
$\alpha_{\mathrm{CE}}\sim 0.2-0.4$, consistent with the range of
$\alpha_{\mathrm{CE}}$ derived for the sample of close WD+BD systems with
accurate parameters (Zorotovic & Schreiber, 2022).

Figure 1: Total age of the system (top), initial mass of the progenitor of the
white dwarf (middle), and orbital period at the onset of the common-envelope
phase (bottom), as a function of the common-envelope efficiency
$\alpha_{\mathrm{CE}}$, for the possible progenitors of the inner white dwarf
in our reconstruction.

##  3 Possible impact of triple dynamics on the evolution

Given that the close inner WD+BD binary is part of a hierarchical triple, it
might in principle be possible that ZLK oscillations had an impact on the
formation of the close binary. If ZLK oscillations were present after the
formation of the inner white dwarf, one cannot a priori exclude that the inner
binary perhaps did not form through common envelope evolution but that instead
ZLK oscillations generated large eccentricities and subsequent tidal decay,
producing the short orbital period of the WD+BD inner binary we observe today.

In what follows, we use the methodology developed by Muñoz & Petrovich (2020)
to evaluate whether the formation of the close WD+BD binary can be understood
by inward migration due to ZLK oscillations coupled with tidal friction after
the formation of the inner white dwarf (henceforth the WD/BD + WD phase).

###  3.1 Constraints for migration during the WD/BD+WD phase

After the formation of the second (inner) white dwarf, migration of the brown
dwarf is achieved if the inner eccentricity during ZLK oscillations is above a
critical value $e_{\mathrm{mig}}$ required for tidal migration, but below
$e_{\mathrm{dis}}$ to avoid tidal disruption. Both $e_{\mathrm{mig}}$ and
$e_{\mathrm{dis}}$ can be obtained analytically as:

$e_{\mathrm{mig}}=1-1.96\left(\frac{k_{\mathrm{2,BD}}\Delta
t_{L}T}{P^{2}_{\mathrm{BD}}}\frac{M_{\mathrm{WD,inner}}}{M_{\mathrm{BD}}}\frac{R^{5}_{\mathrm{BD}}}{a^{5}_{\mathrm{BD}}}\right)^{1/7}$
(1)

and

$e_{\mathrm{dis}}=1-\eta_{\mathrm{dis}}\frac{R_{\mathrm{BD}}}{a_{\mathrm{BD}}}\left(\frac{M_{\mathrm{WD,inner}}}{M_{\mathrm{BD}}}\right)^{1/3}$
(2)

(Eqs. 9 and 10 of Muñoz & Petrovich 2020). Here $M_{\mathrm{BD}}$,
$R_{\mathrm{BD}}$, $k_{\mathrm{2,BD}}$, $a_{\mathrm{BD}}$, and
$P_{\mathrm{BD}}$ are the mass, radius, potential Love number of degree 2,
semi-major axis and orbital period of the brown dwarf, respectively. $\Delta
t_{L}$ is the lag time, $T$ the time interval in which the migration occurs
and $\eta_{\mathrm{dis}}=2.7$ a numerical factor to estimate the minimum
orbital separation at which tidal disruption will occur (Guillochon et al.,
2011).

For the inner white dwarf we assume a progenitor of
$1.07\,\mbox{M${}_{\odot}$}$, which ends up as a white dwarf with
$M_{\mathrm{WD,inner}}\approx 0.52\,\mbox{M${}_{\odot}$}$ after $\approx 9.7$
Gyr according to the SSE stellar evolution code (Hurley et al., 2000). In the
same way, we assume for the outer white dwarf a progenitor of
$1.7\,\mbox{M${}_{\odot}$}$ which formed a $0.6\,\mbox{M${}_{\odot}$}$ white
dwarf after $\approx 2.1$ Gyr. For the brown dwarf mass we use the value
estimated by Rebassa-Mansergas et al. (2022) of $0.07\,\mbox{M${}_{\odot}$}$
while its radius is assumed to be $0.086\,\mbox{R${}_{\odot}$}$ based on the
brown dwarf isochrones222http://perso.ens-
lyon.fr/isabelle.baraffe/COND03_models of Baraffe et al. (2003).
$k_{\mathrm{2,BD}}$ is set to $0.286$ (Heller et al., 2010), $\Delta t_{L}=1$
s (i.e. 10 times the Jupiter lag-time) and $T=0.3$ Gyr assuming a total age of
$10$ Gyr. Figure 2 shows the range of eccentricities
$e_{\mathrm{mig}}<e<e_{\mathrm{dis}}$ that allow for migration as function of
$a_{\mathrm{BD}}$. The migration windows becomes very narrow at
$a_{\mathrm{BD}}\simeq 1$ au, ranging from $0.988\lesssim e\lesssim 0.997$.

During ZLK oscillations, extremely high eccentricities ($e\approx 1$) can be
achieved at the octupolar level of approximation. The importance of this
regime is usually measured by

$\epsilon_{\mathrm{oct}}=\left(\frac{M_{\mathrm{WD,inner}}-M_{\mathrm{BD}}}{M_{\mathrm{WD,inner}}+M_{\mathrm{BD}}}\right)\left(\frac{a_{\mathrm{BD}}}{a_{\mathrm{outer}}}\right)\frac{e_{\mathrm{outer}}}{1-e_{\mathrm{outer}}^{2}},$
(3)

where $a_{\mathrm{outer}}$ and $e_{\mathrm{outer}}$ are the semi-major axis
and eccentricity of the tertiary companion. In general, it is assumed that the
octupole regime plays an important role in the evolution of the system when
$\epsilon_{\mathrm{oct}}\gtrsim 0.001$. Taking a conservative value of $1700$
au for $a_{\mathrm{outer}}$ (based on the projected separation of $\approx
1673$ au between both white dwarfs) and $e_{\mathrm{outer}}=0.5$, the octupole
regime in Gaia 0007$-$1605 becomes important when $a_{\mathrm{BD}}\gtrsim 3$
au.

To estimate the maximum eccentricity $e_{\mathrm{max}}$ attained by the inner
binary in the octuplar regime we use equation (7) of Muñoz & Petrovich (2020),
which is derived from the perpendicular test particle quadrupole (TPQ)
approximation (Liu et al., 2015; Naoz, 2016) and assuming an initial circular
inner orbit:

$\begin{split}0=\frac{9}{8}e^{2}_{\mathrm{max}}-\xi_{\mathrm{GR}}\left(\frac{1}{(1-e^{2}_{\mathrm{max}})^{1/2}}-1\right)\\\
-\frac{\xi_{\mathrm{tide}}}{15}\left(\frac{1+3e^{2}_{\mathrm{max}}+\frac{3}{8}e^{4}_{\mathrm{max}}}{(1-e^{2}_{\mathrm{max}})^{9/2}}-1\right),\end{split}$
(4)

where

$\xi_{\mathrm{GR}}=\frac{3GM^{2}_{\mathrm{WD,inner}}a^{3}_{\mathrm{outer}}(1-e^{2}_{\mathrm{outer}})^{3/2}}{a^{4}_{\mathrm{BD}}c^{2}M_{\mathrm{WD,outer}}}$
(5)

and

$\xi_{\mathrm{tide}}=\frac{15M^{2}_{\mathrm{WD,inner}}a^{3}_{\mathrm{outer}}(1-e^{2}_{\mathrm{outer}})^{3/2}k_{\mathrm{2,BD}}}{a^{8}_{\mathrm{BD}}M_{\mathrm{BD}}M_{\mathrm{WD,outer}}}$
(6)

are terms that represent the strength of general relativistic (GR) precession
and tides relative to the quadrupolar potential of the triple system, with $c$
being the speed of light. Solutions of Equation 4 for the current
configuration of the system are shown in Figure 2. The value of
$e_{\mathrm{max}}$ (magenta line) reaches the migration window when
$a_{\mathrm{BD}}\gtrsim 5.11$ au. The eccentricities for migration correspond
to a small range of rather extreme values ($e_{mig}=0.9984\lesssim
e_{max}\lesssim e_{dis}=0.9998$).

The range of mutual inclinations [$90\degree-\Delta i$, $90\degree+\Delta i$]
in which $e_{\mathrm{max}}$ is attained can be estimated by

$\Delta
i=2.9\degree\left(\frac{\epsilon_{\mathrm{oct}}}{10^{-3}}\right)^{1/2}$ (7)

(equation 6 of Muñoz & Petrovich, 2020). For a brown dwarf located at $5.11$
au prior to migration, $\Delta i\approx 3.6\degree$. If the mutual inclination
is taken from an isotropic distribution (i.e. uniform in $\cos{i}$), then the
probability of Gaia 0007$-$1605 being in the inclination window required for
migration is $\approx 6\%$. It is worth mentioning that the value of $\Delta
i$ can be increased by assuming a larger inner semi-major axis but at the
expense of reducing the eccentricity window for migration (unless the planet
is able to survive the migration when $e_{\mathrm{max}}>e_{\mathrm{dis}}$,
being partially stripped).

If migration due to ZLK oscillations is successful, the final semi-major axis
of the brown dwarf can be approximated as:

$a_{\mathrm{BD,f}}\approx 2a_{\mathrm{BD,0}}(1-e_{\mathrm{max}}),$ (8)

where subscripts f and 0 stand for after and before migration, respectively.
We found that when Equation 8 is evaluated for solutions of equation 4 located
in the migration window, the resulting range of final semi-major axis
($0.0043\lesssim a_{\mathrm{BD,f}}\lesssim 0.0166$ au) is $\approx 1\sigma$
below the one derived using $M_{\mathrm{WD,inner}}=0.52\pm 0.02$,
$M_{\mathrm{BD}}=0.07\pm 0.01$ and the orbital separation calculated from the
orbital period reported by Rebassa-Mansergas et al. (2022, $0.0169\pm
2.8\times 10^{-4}$ au). This result is depicted in the zoomed portion of
Figure 2, in which the eccentricity required (according to equation 8, red
dotted line) for the brown dwarf to reach the observed orbital period is
always below the migration window. This rather small discrepancy can be solved
either by assuming a slightly lower total mass for the white dwarf/brown dwarf
binary (especially a smaller brown dwarf mass than what we assumed can not be
excluded given the weak current observational constraints) or by assuming in
our model a smaller outer semi-major axis so the maximum inner eccentricity
reaches the migration window at smaller inner semi-major axes where solutions
of Equation 8 also are in this window. An example of the latter is depicted by
the green line in Figure 2 which shows $e_{\mathrm{max}}$ when
$a_{\mathrm{outer}}=1115$ au (this value is obtained if we assume
$e_{\mathrm{outer}}=0.5$ and the observed projected separation being the outer
apoapsis distance). In this configuration, $e_{\mathrm{max}}$ reaches the
migration window when $a_{\mathrm{BD}}\gtrsim 3.2$ au and the minimum
eccentricity for migration is $\approx 0.997$ with a similar range of mutual
inclinations.

The results presented above suggest that ZLK oscillations are able to induce
inward migration although two main conditions are required to make this
possible. One needs (1) the inner eccentricity to be at least $\approx 0.997$
which in turn implies that (2) the mutual orbital inclination should be
between $\approx 86\degree$ and $\approx 94\degree$.

Figure 2: Eccentricity window (filled area) in which migration would occur
during the WD/BD+WD phase as function of the semi-major axis of the brown
dwarf. The boundaries of this window are set by Equations 1 (dashed line) and
2 (dash-dotted line). The maximum eccentricity attained by the brown dwarf
under the perpendicular TPQ approximation (i.e. solution of Equation 4) as
function of the brown dwarf semi-major axis and taking
$a_{\mathrm{WD,outer}}=1115,1700$ and $2500$ au (green, magenta and brown
lines respectively). The eccentricity required for the brown dwarf to reach
its current orbital period (semi-major axis) after successful migration is
given by the red dotted line which has been calculated using Equation 8.

###  3.2 Validity of the analytical model used

The analytical model used by Muñoz & Petrovich (2020), developed to find the
parameter space in which the tight orbit of the planet candidate in WD 1856
can be explained by tidal migration due to ZLK oscillations, concluded that
such configuration is only reproduced if the initial (i.e. during the main-
sequence stage of the host and tertiary) semi-major axis of the planet is
located within a narrow range of values. This result, however, is in
disagreement with the outcome of the numerical simulations performed by
Stephan et al. (2021), who derived a much wider range of initial semi-major
axis. As discussed by Stephan et al. (2021), this discrepancy can be explained
considering two differences in the modeling: (1) The analytical condition
imposed by Muñoz & Petrovich (2020) for the planet to survive the evolution of
its host are much more restrictive than the numerical simulations performed by
Stephan et al. (2021), (2) unlike Muñoz & Petrovich (2020), Stephan et al.
(2021) consider the scenario in which the planet crosses the Roche limit (i.e.
the planet reaches eccentricities above $e_{\mathrm{dis}}$) and is partially
stripped by tidal forces but survives the migration process.

Because of these limitations of the analytical prescriptions from Muñoz &
Petrovich (2020), in our analysis we simply assumed that the brown dwarf
survived the evolution of its host and emphasize only the minimum eccentricity
required to achieve migration. In addition, unlike Muñoz & Petrovich (2020) we
treated the outer semi-major axis as a free parameter. In this regard it is
important to recall that the outer semi-major axis used in section 3.1 serves
as an example and does not exclude configurations with smaller (larger)
values, which will move the maximum eccentricity attained by the brown dwarf
towards smaller (larger) inner semi-major axes. For instance, taking
$a_{\mathrm{outer}}=2500$ au the minimum eccentricity for migration is
$\approx 0.9989$ at $a_{\mathrm{BD}}\approx 7$ au as shown by the brown line
in Figure 2.

This result is consistent with the findings of Stephan et al. (2021) on WD
1856 which has a very similar configuration. The range of initial semi-major
axes Stephan et al. (2021) derived for the planet orbiting WD 1856 peaks close
to 100 au. At such distances, $\xi_{\mathrm{tide}}$ and $\xi_{\mathrm{GR}}$
become negligible and the maximum eccentricity only depends on the mutual
inclination $i$ between the inner and outer orbits
($e_{\mathrm{max}}=(1-5\cos^{2}{(i)}/3)^{1/2}$). According to Equation 8, and
neglecting any orbital expansion due to stellar evolution mass loss333This
assumption provides a lower limit on the true eccentricity required for
migration. If the orbit expands adiabatically (i.e. the orbital separation
increases by a factor $M_{\mathrm{tot,0}}/M_{\mathrm{tot,f}}$, where
$M_{\mathrm{tot}}$ is the total mass of the inner binary and subscripts f and
0 stand for after and before mass loss) the brown dwarf/planet semi-major axis
would increase by a factor of $\approx 1.9$ after the formation of the host
white dwarf., a brown dwarf/planet located at 100 au from its host star/white
dwarf would require an eccentricity of $\approx 0.9999$ to migrate to 0.02 au,
which in turn would require $i\approx 89.4\degree$.

###  3.3 Migration scenario due to ZLK oscillations coupled with common-
envelope evolution

Figure 3: Minimum eccentricity required for RLOF (Equation 9) as function of
the mass of the inner white dwarf progenitor during the AGB from the SSE
model. For a brown dwarf semi-major axis equal to $5$ au ($10$ au) the minimum
eccentricity is $\approx 0.62$ ($\approx 0.81$) when the host star has a mass
of $\approx 0.67$ M⊙ on the thermally pulsating AGB.

Even if ZLK oscillations did not produce eccentricities high enough for tidal
migration, lower values of the eccentricity could still have impacted the
evolution and perhaps triggered common envelope evolution, similar to what has
been proposed for the transiting planet around WD 1856+534 (e.g. Lagos et al.,
2021; Trani et al., 2022).

In order to estimate the minimum eccentricity required to start Roche lobe
overflow (RLOF) we use the Roche lobe radius approximation for eccentric
orbits given by Sepinsky et al. (2007, their Equation 45) evaluated at the
periapsis of the binary and taking the Roche radius equal to the stellar
radius $R_{\star}$ at the onset of mass transfer:

$e_{\mathrm{min,RLOF}}=1-\frac{R_{\star}}{a_{\mathrm{BD}}}\frac{0.6q^{2/3}+\ln{(1+q^{1/3})}}{0.49q^{2/3}}.$
(9)

Here $q=M_{\mathrm{\star}}/M_{\mathrm{BD}}$ is the mass ratio between the host
star and the brown dwarf. Both $R_{\star}$ and $M_{\mathrm{\star}}$ are
obtained from SSE. Figure 3 shows that for a separation of 5 (10) au, the
eccentricity required to trigger mass transfer is
$e_{\mathrm{min,RLOF}}\approx 0.62$ ($\approx 0.81$) when the host star has a
mass of $\approx 0.67$ M⊙ during the AGB phase.

To provide one exemplary full simulation of the outlined triple star
evolution, we used the Multiple Stellar
Evolution444https://github.com/hamers/mse. MSE is based on the SSE (Hurley et
al., 2000) and binary stellar evolution (BSE, Hurley et al., 2002) codes.
(MSE, Hamers et al., 2021) code version 0.86, which allows to calculate the
secular orbital evolution of the brown dwarf including the effects of stellar
and tidal evolution, general relativity, N-body dynamics and binary
interactions. Table 2 summarizes the initial parameters used in our
simulation. After $\approx 2$ Gyr the tertiary star evolves into a white dwarf
and its orbit expands, increasing the timescale of the ZLK oscillations but
keeping approximately the same maximum eccentricity. When the host star
evolves through the AGB, the inner eccentricity is still high enough
($e\approx 0.73$) to trigger RLOF. Although this simulation assumes an outer
semi-major axis of $\approx 700$ au (which according to the adiabatic mass
loss model will increase to the observed projected separation of $\sim 1700$
au ), we verified with additional simulations that for an initial outer semi-
major axis up to $1200$ au ZLK oscillations with tidal friction are still able
to lead to RLOF during the AGB. It therefore appears to be possible that the
ZLK mechanism played a role in the evolution of the triple system prior to
common envelope evolution.

Figure 4: Orbital evolution of Gaia 0007$-$1605 according to the initial
orbital configuration given in Table 2. The simulation is stopped when RLOF
starts in the inner binary (red circle). Top panel: Evolution of the brown
dwarf eccentricity. While the host and tertiary stars are on the main sequence
($\lesssim 2$ Gyr) ZLK oscillations produce a maximum inner eccentricity of
$\approx 0.85$. After $\approx 2$ Gyr the tertiary star evolves into a white
dwarf, its semi-major axis increases and therefore also the timescale of the
ZLK oscillations. Bottom panel: Evolution of the tertiary star semi-major axis
and periapsis (solid and dashed magenta lines respectively), brown dwarf semi-
major axis (green line), stellar radius and Roche radius at periapsis of the
host star (red and black lines respectively). At $\approx 9.7$ Gyr the
eccentricity of the inner binary is high enough ($\approx 0.73$) to trigger
RLOF (denoted by the red circle) when the semi-major axis of the brown dwarf
and the tertiary star are $\approx 1.5$ au (orbital period of $\approx 720$ d)
and $\approx 1400$ au respectively. Table 2: Configuration assumed for Gaia
0007$-$1605 at the beginning of the MSE simulation. Stellar radii and apsidal
motion constants ($k_{\mathrm{AM}}$) are set internally by the code. For the
brown dwarf we use $R_{\mathrm{BD}}=0.086$ R⊙(expressed in au) and
$k_{\mathrm{AM}}=0.143$ (i.e. half of the Love number of degree 2). The number
of output steps is set to 20 000. The mutual inclination between the inner and
outer orbits is $65.8\degree$.

Parameter | Host star | Tertiary star | Brown dwarf
---|---|---|---
Mass [M⊙] | 1.07 | 1.7 | 0.07
Radius [au] | Default | Default | $4\times 10^{-3}$
Metallicity | Solar | Solar | Solar
Apsidal motion constant | 0.19 | 0.19 | 0.143
Orbital parameter | Inner orbit | Outer orbit |
Eccentricity | 0 | 0.45 |
Semi-major axis [au] | 4.5 | 720 |
Inclination [rad] | 0.001 | 1.15 |
Argument of pericentre [rad] | 0.01 | 0.5 |
Longitude of ascending node [rad] | 0.5 | 0.5 |

##  4 Conclusion

We have studied in detail the evolutionary history of the WD+BD binary Gaia
0007$-$1605 which is the inner binary of a hierarchical triple star system
with the tertiary being a white dwarf that is $\approx 8$ Gyr older than the
inner one. We found that assuming the planet survived the evolution of its
host star into a white dwarf, ZLK oscillations alone can explain the
configuration we observe today if the inner eccentricity reached values close
or above $\approx 0.997$, which in turn implies that the inner and outer
orbits were (or are) close to being perpendicular to each other.

By reconstructing the close orbit of the brown dwarf we observe today through
common envelope evolution we found that no energy in addition to orbital
energy is required to understand the currently observed period. Our findings
further support the conclusions recently drawn by Zorotovic & Schreiber (2022)
that, in contrast to previous findings (e.g. De Marco et al., 2011), common
envelope evolution with substellar companions does not require additional
energy sources to play a role. This also further illustrates that in terms of
the energy budget of common envelope evolution, the transiting planet around
WD 1856+534 remains an outlier as additional energy is required to reproduce
the currently observed period with common envelope evolution. The idea that an
additional planet might have contributed to the ejection of the envelope (e.g.
Bear & Soker, 2011), as suggested first by Lagos et al. (2021) and later
investigated in more detail by Chamandy et al. (2021), offers an elegant
solution for this particular system (although it will be impossible to verify
this hypothesis).

However, when considering constraints on common envelope evolution from WD+BD
binaries it is important to keep in mind that Gaia 0007$-$1605 and WD 1856+534
might not be post common envelope systems, i.e. both common envelope evolution
and/or ZLK oscillations plus tidal migration appear as potential scenarios to
have produced the currently observed configuration.

## Acknowledgements

We thank the referee for helpful comments and suggestions. This project has
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (Grant agreement No.
101020057). MRS and MZ acknowledge support from FONDECYT (grant 1221059). MRS
also acknowledges support by ANID, – Millennium Science Initiative Program –
NCN19_171. This research was supported in part by the National Science
Foundation under Grant No. NSF PHY-1748958 (KITP).

## Data Availability

The data and numerical tools used in this article can be obtained upon request
to the corresponding author and after agreeing to the terms of use.

## References

  * Alonso et al. (2021) Alonso R., et al., 2021, A&A, 649, A131
  * Baraffe et al. (2003) Baraffe I., Chabrier G., Barman T. S., Allard F., Hauschildt P. H., 2003, A&A, 402, 701
  * Bear & Soker (2011) Bear E., Soker N., 2011, MNRAS, 411, 1792
  * Chamandy et al. (2021) Chamandy L., Blackman E. G., Nordhaus J., Wilson E., 2021, MNRAS, 502, L110
  * De Marco et al. (2011) De Marco O., Passy J.-C., Moe M., Herwig F., Mac Low M.-M., Paxton B., 2011, MNRAS, 411, 2277
  * El-Badry et al. (2021) El-Badry K., Rix H.-W., Heintz T. M., 2021, MNRAS, 506, 2269
  * Gentile Fusillo et al. (2021) Gentile Fusillo N. P., et al., 2021, MNRAS, 508, 3877
  * Guillochon et al. (2011) Guillochon J., Ramirez-Ruiz E., Lin D., 2011, ApJ, 732, 74
  * Hamers et al. (2021) Hamers A. S., Rantala A., Neunteufel P., Preece H., Vynatheya P., 2021, MNRAS, 502, 4479
  * Heller et al. (2010) Heller R., Jackson B., Barnes R., Greenberg R., Homeier D., 2010, A&A, 514, A22
  * Hernandez et al. (2021) Hernandez M. S., et al., 2021, MNRAS, 501, 1677
  * Hernandez et al. (2022a) Hernandez M. S., et al., 2022a, MNRAS, 512, 1843
  * Hernandez et al. (2022b) Hernandez M. S., et al., 2022b, MNRAS, 517, 2867
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Ivanova et al. (2013) Ivanova N., et al., 2013, A&ARv, 21, 59
  * Kozai (1962) Kozai Y., 1962, AJ, 67, 591
  * Lagos et al. (2021) Lagos F., Schreiber M. R., Zorotovic M., Gänsicke B. T., Ronco M. P., Hamers A. S., 2021, MNRAS, 501, 676
  * Lagos et al. (2022) Lagos F., Schreiber M. R., Parsons S. G., Toloza O., Gänsicke B. T., Hernandez M. S., Schmidtobreick L., Belloni D., 2022, MNRAS, 512, 2625
  * Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
  * Liu et al. (2015) Liu B., Muñoz D. J., Lai D., 2015, MNRAS, 447, 747
  * Masuda et al. (2019) Masuda K., Kawahara H., Latham D. W., Bieryla A., Kunitomo M., MacLeod M., Aoki W., 2019, ApJ, 881, L3
  * Merlov et al. (2021) Merlov A., Bear E., Soker N., 2021, ApJ, 915, L34
  * Muñoz & Petrovich (2020) Muñoz D. J., Petrovich C., 2020, ApJ, 904, L3
  * Naoz (2016) Naoz S., 2016, ARA&A, 54, 441
  * Naoz et al. (2011) Naoz S., Farr W. M., Lithwick Y., Rasio F. A., Teyssandier J., 2011, Nature, 473, 187
  * Nelemans et al. (2000) Nelemans G., Verbunt F., Yungelson L. R., Portegies Zwart S. F., 2000, A&A, 360, 1011
  * O’Connor et al. (2021) O’Connor C. E., Liu B., Lai D., 2021, MNRAS, 501, 507
  * Paczynski (1976) Paczynski B., 1976, in Eggleton P., Mitton S., Whelan J., eds, IAU Symposium Vol. 73, Structure and Evolution of Close Binary Systems. p. 75
  * Perets & Fabrycky (2009) Perets H. B., Fabrycky D. C., 2009, ApJ, 697, 1048
  * Petrovich (2015) Petrovich C., 2015, ApJ, 799, 27
  * Rebassa-Mansergas et al. (2021) Rebassa-Mansergas A., et al., 2021, MNRAS, 505, 3165
  * Rebassa-Mansergas et al. (2022) Rebassa-Mansergas A., Xu S., Raddi R., Pala A. F., Solano E., Torres S., Jiménez-Esteban F., Cruz P., 2022, ApJ, 927, L31
  * Schreiber & Gänsicke (2003) Schreiber M. R., Gänsicke B. T., 2003, A&A, 406, 305
  * Sepinsky et al. (2007) Sepinsky J. F., Willems B., Kalogera V., 2007, ApJ, 660, 1624
  * Stephan et al. (2021) Stephan A. P., Naoz S., Gaudi B. S., 2021, ApJ, 922, 4
  * Thompson (2011) Thompson T. A., 2011, ApJ, 741, 82
  * Trani et al. (2022) Trani A. A., Rieder S., Tanikawa A., Iorio G., Martini R., Karelin G., Glanz H., Portegies Zwart S., 2022, arXiv e-prints, p. arXiv:2205.13537
  * Vanderburg et al. (2020) Vanderburg A., et al., 2020, arXiv e-prints, p. arXiv:2009.07282
  * Von Zeipel (1910) Von Zeipel H., 1910, Astronomische Nachrichten, 183, 345
  * Webbink (1984) Webbink R. F., 1984, ApJ, 277, 355
  * Webbink (2008) Webbink R. F., 2008, Common Envelope Evolution Redux. In: Short-Period Binary Stars: Observations, Analyses, and Results. Astrophysics and Space Science Library, vol 352. Springer, Dordrecht, p. 233, doi:10.1007/978-1-4020-6544-6_13
  * Woods et al. (2012) Woods T. E., Ivanova N., van der Sluys M. V., Chaichenets S., 2012, ApJ, 744, 12
  * Wu & Murray (2003) Wu Y., Murray N., 2003, ApJ, 589, 605
  * Xu et al. (2021) Xu S., et al., 2021, AJ, 162, 296
  * Zorotovic & Schreiber (2022) Zorotovic M., Schreiber M., 2022, MNRAS, 513, 3587
  * Zorotovic et al. (2010) Zorotovic M., Schreiber M. R., Gänsicke B. T., Nebot Gómez-Morán A., 2010, A&A, 520, A86


A quantum annealing approach to minimum distance problem

Refat Ismail
Department of Physics and Astronomy, University of Kentucky, Lexington, KY, 40503 USA

Ashish Kakkar

Anatoly Dymarsky
Department of Physics and Astronomy, University of Kentucky, Lexington, KY, 40503 USA

Quantum error-correcting codes (QECCs) is at the heart of fault-tolerant quantum computing.
As the size of quantum platforms is expected to grow, one of the open questions is to design new optimal  codes of ever-increasing size. A related challenge is to “certify” the quality of a given code by evaluating its minimum distance, a quantity characterizing code's capacity to preserve quantum information. This problem is known to be NP-hard. Here we propose to harness the power of contemporary quantum platforms to address this question, and in this way  help design quantum platforms of the future.
Namely, we introduce an approach to compute the minimum distance of quantum stabilizer codes by reformulating the problem as a Quadratic Unconstrained Binary Optimization (QUBO) problem and leveraging established QUBO algorithms and heuristics as well as quantum annealing (QA) to address the latter. The reformulation as a QUBO introduces only a logarithmic multiplicative overhead in the required number of variables. We demonstrate practical viability of our method by comparing the performance of purely classical algorithms with the D-Wave Advantage 4.1 quantum annealer as well as hybrid quantum-classical algorithm Qbsolv. We found that the hybrid approach demonstrates competitive performance, on par with the best available classical algorithms to solve QUBO. In a practical sense the QUBO-based approach is currently lagging behind the best deterministic minimal distance algorithms, however this advantage may disappear as the size of the platforms grows.

§ INTRODUCTION

The field of quantum error correction has been the focus of intense research in recent decades due to its role in the development of fault-tolerant quantum computing. Quantum fault tolerance offers a trade-off between efficiency in storing and manipulating information and robust protection against errors <cit.>. One of the central theoretical goals is to design QECCs that minimize the required resources while enhancing the protective properties.

Quantum algorithms operate on an abstract Hilbert space, known as the logical Hilbert space. When employing a QECC, we restrict the logical Hilbert space to a specific subspace of the full computational (also referred to as physical) Hilbert space. Errors that send the computation outside the code subspace could be detected and corrected. In contrast, errors that map the code subspace into itself are undetectable. The code's minimum (or Hamming) distance $d$ is defined as the smallest weight of an undetectable error, i.e. how many physical qubits this  error affects <cit.>.

A QECC can be characterized by three parameters: the number of physical qubits $n$, the number of logical states $K$ it protects and its minimum distance $d$. Corresponding codes are conventionally denoted as $((n, K, d))$. Constructing QECCs of ever-increasing $n$, with a fixed ratio $\ln K/n$ and the best (largest) $d/n$ is an open problem. Our emphasis will be on stabilizer codes, the most common class of quantum codes introduced in <cit.> and widely considered the best strategy to achieve fault tolerance in the near term.
A binary stabilizer code, i.e. over qubits, protects the $K=2^k$ states, where $k$ denotes the number of logical qubits and is denoted by $[[n,k,d]]_2$.

An optimal  code is the one with the highest known value of $d$ for given $n$ and $k$. In practice to certify that code is optimal or close to optimal is a difficult problem.
Depending on how the code is defined, determining $d$ could be a computationally challenging task. In fact, finding $d$ for a general (binary) stabilizer code specified by its generating matrix is known to be NP-hard. Even approximating the minimum distance, with a guaranteed precision specified by a given multiplicative or additive factor, is also NP-hard <cit.>.  There exist multiple approaches for computing the minimum distance of quantum codes.  One of the best deterministic algorithms was developed by G. White and M. Grassl and has been incorporated into MAGMA <cit.>. Multiple additional techniques have been also tailored for specific families of quantum codes, such as  Low-Density-Parity-Check (qLDPC) codes <cit.>.

In this paper, we introduce a new computational technique to evaluate or bound from above the minimum distance of binary quantum stabilizer codes, by reformulating the question as a Quadratic Unconstrained Binary Optimization (QUBO) problem. A QUBO is a general framework to represent combinatorial optimization problem as a problem to minimize quadratic functions of binary variables <cit.>. To reformulate the minimum distance problem as a QUBO, one needs to introduce auxiliary variables, e.g. to implement mod $2$ arithmetic. As a result, our method involves a multiplicative logarithmic overhead in terms of the total number of variables.

Reformulating the problem in terms of QUBO introduces a promising pathway to harness power of quantum computational schemes such  as Adiabatic Quantum Computing (AQC) and Quantum Annealing (QA), alongside a spectrum of established classical algorithms and heuristics designed specifically for QUBO problems. There is an ongoing effort to reformulate many NP-hard problems of theoretical and practical value within a unified framework of QUBO, with such problems as SAT, Knapsack, Maximum Cut, and Graph Coloring already being worked out  in the literature <cit.>. These problems arise in different fields, from Quantum Information to Machine Learning to RNA folding <cit.>. Our work can be seen as contributing to  this  effort.

To investigate the efficiency of our approach, implemented on the contemporary quantum simulators, we focus on quantum annealing (QA) as a method to solve QUBO. QA is a quantum heuristic approach, inspired by the principles of Adiabatic Quantum Computing (AQC) <cit.>. QA provides a rigorous upper bound for optimization problems, albeit without a guarantee on how close the obtained result is to the optimal value.
To investigate the viability of QA in our case, we tested our method using both the D-Wave Advantage 4.1 quantum annealer and Qbsolv, a hybrid quantum-classical algorithm developed by D-Wave and running on their platform <cit.>. As a benchmark, we compared the performance of quantum and quantum-classical methods with simulated annealing (SA), a well-established classical heuristic approach <cit.>. Whenever available, we compared the results of quantum, classical, or hybrid optimization with the exact results.

Our study has shown that the hybrid quantum/classical approach is already demonstrating performance at the level of purely classical optimization scheme, possibly opening a new avenue for near-future practical applications of quantum annealing.

The paper is organized as follows. In Section II, we introduce the basics of stabilizer codes. We formulate the main algorithm to evaluate $d$ as a QUBO problem in Section III. Section IV provides a brief overview of AQC and QA. Section V details the comparison of purely classical, quantum, and hybrid approaches. An open-source implementation of our algorithm is available on GitHub <cit.>. We conclude with a discussion in section VI.

§ THEORETICAL PRELIMINARIES

§.§ Stabilizer Codes: a brief overview

One of the most frequently used  class of QECCs is the stabilizer codes. If S is an Abelian subgroup of the Pauli group $\mathcal{P}_n$[The $n$ qubit Pauli group, $\mathcal{P}_n$,  consists of all n-fold tensor products of  the Pauli operators $I, X, Y, Z$, with an overall phase $i^k$ with an integer $k$.] with $(n-k)$ generators and $-I \notin S$, then the simultaneous +1 eigenspace, $V_S$, of all elements of $S$ is an $[[n,k,d]]$ quantum stabilizer code, where $n$, $k$, and $d$ denote the number of physical qubits, logical qubits, and the code minimum distance, correspondingly. Note that the two conditions on the subgroup, being Abelian and $-I \notin S$, are to ensure that the code space is non-trivial.

To specify the stabilizer code, it is sufficient to provide $(n-k)$ generators of the stabilizer group $S$. The stabilizer generators contain essential information about the error-correcting capabilities of the code, eliminating the need to construct the code space $V_S$ explicitly. Stabilizer codes are also known as additive quantum codes due to their one-to-one correspondence with trace hermitian self-orthogonal additive classical codes over GF(4) <cit.>.

Pauli group elements can be represented as binary vectors, allowing for an efficient representation of the stabilizer group. Namely, up to a phase, any Pauli group element can be expressed as a product of X and Z operators, each raised to a power of either 0 or 1,

\begin{align}
g \propto \left( X_1^{\alpha_1} X_2^{\alpha_2}\ldots X_n^{\alpha_n} \right) \cdot \left( Z_1^{\beta_1} Z_2^{\beta_2}\ldots Z_n^{\beta_n} \right), \quad g \in \mathcal{P}_n
\label{eq:stabilizerfromvec}
\end{align}

Here $\alpha_i$ and $\beta_i$ are binary variables.  Combining  them together, we can represent a Pauli element using a binary vector $\boldsymbol{c}=(\boldsymbol{\alpha}|\boldsymbol{\beta})$[We use bold symbols to represent vectors.]. Multiplying two Paulis then corresponds to adding their respective binary vectors mod 2.

In this representation, the stabilizer group generators can be specified by a binary $(n-k) \times 2n$ matrix, where each row is vector representation of one of the $(n-k)$ generators. This matrix H is called the stabilizer matrix (also known as the parity check matrix). We can generate  binary representations of all elements of the stabilizer group as linear combinations (mod 2) of the generators of the group,

\begin{equation}
\boldsymbol{c} = H^T\boldsymbol{x} \ \text{(mod 2)}, \quad \boldsymbol{x} \in \mathbb{Z}^{n-k}_2,  \quad \boldsymbol{c} \in \mathbb{Z}^{2n}_2,\quad \mathbb{Z}_2\equiv\{0,1\},
\end{equation}

where $\boldsymbol{c}$ is the binary representation of a given stabilizer group element, sometimes called a  “codeword.”

The normalizer of S, denoted $N_{\mathcal{P}_n}(S)$, consists of all elements of Pauli group that map the stabilizer group to itself under conjugation, $N_{\mathcal{P}_n}(S) = \{g \in \mathcal{P}_n \ | \ gSg^{-1}= S \}$. It is clear that $S \subseteq N_{\mathcal{P}_n}(S)$.  The operators present in the normalizer but not in stabilizer represent logical Pauli operators on the $k$ encoded logical qubits. These operators map states from the code subspace into each other, which makes them undetectable errors.

The generators of the normalizer group, $N_{\mathcal{P}_n}(S)$, can be represented as columns of the normalizer matrix G (also known as the generator matrix), which is a $(2n) \times (n+k)$ matrix of full rank that satisfies

\begin{equation}
\label{HG}
H\, \Lambda \, G = 0 \ (\text{mod 2}), \quad \Lambda :=
\begin{pmatrix}
0 & I_n\\I_n & 0
\end{pmatrix},
\end{equation}

with $I_n$ being the $n \times n$ identity matrix. The full rank of G ensures that all generators of the normalizer are linearly independent, while the condition (<ref>) ensures that all normalizer elements commute with the stabilizer group. Matrix $G$ also satisfies   ${\rm dim}\ker\left(G^T \Lambda G \right) = n-k$.

§.§ Minimum distance

The weight of a Pauli operator – an element from the Pauli group – is the number of qubits on which it acts non-trivially. The weight of a Pauli operator can be calculated using the binary representation $\boldsymbol{c} = (\boldsymbol{\alpha} | \boldsymbol{\beta}) $ as follows,

\begin{equation}
\begin{split}
w(\boldsymbol{c}) = \boldsymbol{\alpha} \cdot \boldsymbol{\alpha} + \boldsymbol{\beta} \cdot \boldsymbol{\beta} - \boldsymbol{\alpha} \cdot \boldsymbol{\beta}.
\end{split}
\label{pauli_weight}
\end{equation}

The minimum distance of a $(k \neq 0)$ stabilizer code is the smallest weight of an operator in the set of all logical Pauli operators, i.e., the set of elements in the normalizer $N_{\mathcal{P}_n}(S)$ but not in the stabilizer $S$. The minimum code distance is also called the minimal Hamming distance.

The code minimum distance $d$ can then be expressed in terms of the normalizer matrix, $G$,

\begin{equation}
d = \min_{\boldsymbol{c}} w(\boldsymbol{c}), \quad \boldsymbol{c} = G\, \boldsymbol{x} \ \text{(mod 2)}, \quad \boldsymbol{x} \in \mathbb{Z}_2^{n+k}, \quad \boldsymbol{c} \notin \text{ker}(G^T\Lambda).
\label{distance_non_dual}
\end{equation}

Here the constraint $\boldsymbol{c} \notin \text{ker}(G^T\Lambda)$ is equivalent to $\boldsymbol{c} \notin \text{Im}(H^T)$ which ensures that we minimize only over the set of logical Pauli operations. A stabilizer code with minimum distance $d$ can correct all errors of weight $\lfloor \frac{d-1}{2} \rfloor$ or less, and can also correct all erasure errors[Erasure errors are errors in which the qubits affected by the error are known.] of weight $d-1$ or less <cit.>.

§.§ Self-dual stabilizer codes

Self-dual codes represent a special class of stabilizer codes with $k=0$. In this case the code space $V_S$  is one-dimensional, often referred to as a stabilizer state. These codes cannot store or be used to manipulate quantum information, instead, they are protocols for error detection.  By initializing a quantum processor in a stabilizer state and evolving in time such that it should return to the original state,  we can detect and identify any errors that may have accumulated during this evolution. For self-dual codes, the minimum distance of the code is the smallest weight among the weights of all non-trivial elements of the stabilizer group,

\begin{equation}
d = \min_{\boldsymbol{c}} w(\boldsymbol{c}), \quad \boldsymbol{c} = G\boldsymbol{x} \ \text{(mod 2)}, \quad \boldsymbol{x} \in \mathbb{Z}_2^{n}, \quad \boldsymbol{c} \neq 0^{2n}.
\label{distance_dual}
\end{equation}

For self-dual codes, we can chose $G = H^T$. An $[[n,0,d]]$ code can detect all errors affecting up to $d-1$ qubits, and can identify all errors of weight $\lfloor \frac{d-1}{2} \rfloor$ or less <cit.>.

§ ALGORITHM

In this section, we reformulate the problem of finding the minimum distance of a stabilizer code as a QUBO
(quadratic unconstrained binary optimization)  problem.

§.§ Quadratic binary cost function

Any generator matrix can be written as $G^T =\begin{pmatrix}    A & B \end{pmatrix} $, where $A$ and $B$ are $n \times (n+k) $ binary matrices. Using the definition of the Pauli operator weight (<ref>), the code minimum distance (<ref>) can be written as

\begin{eqnarray}
\nonumber
d &=& \min_{\boldsymbol{x} \notin \mathcal{C} } \  \sum_{i=1}^{n} \left(  \alpha_i^2 +\beta_i^2 -\alpha_i \beta_i\right),\ \boldsymbol{x}\in \mathbb{Z}_2^{n+k},\, \,  \label{min_problem}\\
\mathcal{C} &=&  \left\{
\begin{array}{rc}
\{0^{n}\},    &\quad k=0,\\
\{\boldsymbol{x}\, |\,  G^T\Lambda G\,\boldsymbol{x} = \boldsymbol{0}\},   & \quad k\neq 0.
\end{array}
\right.
\label{optimization_prob}
\end{eqnarray}

Here, $\boldsymbol{\alpha} = (A\, \boldsymbol{x} \,\,\,\text{mod 2} )$ and $\boldsymbol{\beta} =(B\, \boldsymbol{x} \,\,\,\text{mod 2} )$. The minimization problem  above  is a constrained optimization problem with a quadratic cost function in terms of binary variables $\{x_i\}$. It closely resembles a QUBO problem <cit.>, except for the presence of constraints and mod 2 sums in the definition of $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$. The subsequent sections explain how to deal with each of these  differences.

§.§ Reformulation of the mod 2  summation

Here, we describe a simple method to trade modular sums into conventional ones by introducing auxiliary integer variables, which subsequently can be parameterized by several binary variables. This leads to a quadratic binary optimization problem that involves only conventional sums.

To eliminate modular sums, we introduce two vectors,

\begin{equation}
\boldsymbol{a} := A\, \boldsymbol{x}, \qquad \boldsymbol{b}  := B\, \boldsymbol{x}.
\label{define_variables}
\end{equation}

Here matrix multiplication assumes conventional summation. Next, we introduce parity $\mathbb{P}$,

\begin{equation}
\label{P}
\mathbb{P}(m)=m\, {\rm mod}\, 2.
\end{equation}

With this definition, we readily recognize $\alpha_i = \mathbb{P}(a_i)$ and $\beta_i = \mathbb{P}(b_i)$. It can be checked straightforwardly that

\begin{equation}
\alpha_i^2 + \beta_i^2 - \alpha_i\beta_i = \frac{1}{2} \left( \mathbb{P}(a_i)  + \mathbb{P}(b_i) + \mathbb{P}(a_i + b_i) \right).
\label{parity_identity}
\end{equation}

To implement (<ref>) we introduce  minimization over an auxiliary integer variable $t$,

\begin{equation}
\mathbb{P}(m) = \min_{t}\,  (m -2t)^2, \quad 0 \leq t \leq \lfloor\frac{\max{m}}{2} \rfloor.
\end{equation}

The resulting cost function becomes

\begin{equation}
{\sf d} = \min_{\boldsymbol {t,u,v}}  \frac{1}{2}\left( \sum_{i=1}^n (a_i - 2t_i)^2 +(b_i - 2u_i)^2 + (a_i +b_i - 2v_i)^2 \right),
\end{equation}

where $\boldsymbol{t,u,v}$ are vectors of auxiliary integer variables. To express integer variables in terms of binary variables, we use their binary representation, i.e.

\begin{equation}
t_i = \sum^{s}_{j=1} 2^{j-1} t_{ij}, \quad t_{ij} \in \mathbb{Z}_2.
\end{equation}

Here $s = \lceil\ln{\left(\frac{\max{a_i}}{2} \right)} \rceil$ is the number of
auxiliary binary variables. This number  is bounded by $\ln{\left(\frac{n+k}{2}\right)}$. Putting it all together, we arrive at

\begin{equation}
d =  \min_{\boldsymbol{x}\notin \mathcal{C}} \min_{\boldsymbol{t,u,v}} \left( \begin{array}{cc}\boldsymbol{x}^T\left( A^TA + B^TB + B^TA\right)\boldsymbol{x} \\[12pt]  - \sum\limits_{i=1}^n \sum\limits_{j=1}^{n+k} \frac{x_j}{2}\left( \sum\limits_{l=1}^{s} 2^l\left( A_{ij} t_{il} +B_{ij} u_{il}  \right)+ \sum\limits_{l=1}^{2s}\left(A_{ij}+ B_{ij} \right) 2^{l} v_{il} \right) \\[14pt]
+ \sum\limits_{i=1}^n\sum\limits_{j,l=1}^{s} 2^{j+l-1}\left(t_{ij}t_{il} + u_{ij}u_{il} \right) + \sum\limits_{i=1}^n\sum\limits_{j,l=1}^{2s} 2^{j+l-1}v_{ij}v_{il}  \end{array} \right),
\end{equation}

where $\boldsymbol{t,u,v}$ now represent matrices of  binary variables $t_{ij},u_{ij},v_{ij}$. The cost function above can be expressed in terms of a binary symmetric matrix

\begin{equation}
\label{ddef}
{\sf d} = %\min_{\boldsymbol{z}}
\boldsymbol{z}^TQ\,\boldsymbol{z},
\end{equation}

where $\boldsymbol{z}$ is a vector of the primary and auxiliary binary variables

\begin{equation}
\boldsymbol{z}= \left( x_1,x_2, \dots ,x_{n+k}, \ t_{1,1},t_{2,1}, \dots  ,t_{n,1}, \dots t_{1,s},\dots ,t_{n,s},u_{1,1},\dots ,u_{n,s},v_{1,1},\dots ,v_{n,2s} \right),
\end{equation}

\begin{equation}
Q :=
\begin{pmatrix}   A^TA + B^TB + B^TA & -(\boldsymbol{w_s} \otimes A)^T &  -(\boldsymbol{w_s} \otimes B)^T & -(\boldsymbol{w_{2s}} \otimes (A+B))^T \\ -\boldsymbol{w_s} \otimes A & 8 \boldsymbol{w_s}\boldsymbol{w_s}^T \otimes I_n  & 0 &0\\
-\boldsymbol{w_s} \otimes B & 0 & 8\boldsymbol{w_s}\boldsymbol{w_s}^T \otimes I_n & 0 \\-
\boldsymbol{w_{2s}} \otimes (A+B) & 0 &0 & 8\boldsymbol{w_{2s}}\boldsymbol{w_{2s}}^T \otimes I_n
\end{pmatrix}.
\label{matrix_eq}
\end{equation}

Here we introduced the column vector  $\boldsymbol{w_s}:= \frac{1}{2}(1,2,4,\dots,2^{s-1})$.  The total number of auxillary binary variables is $4 n s$, with $s \leq \ln{\frac{n+k}{2}} <\ln{n}$. Thus, the total number of binary variables, both primary and auxiliary, grows asymptotically as $O \left(n  \ln{n/n_0}   \right)$, for some constant $n_0$. Notably, the multiplicative overhead in the number of variables increases only logarithmically with $n$.
Finding the minimum distance of the corresponding stabilizer code is equivalent to finding the minimum of (<ref>) subject to constraints.

§.§.§ Self dual codes

For self-dual codes, we can leverage code equivalence (discussed in the next section) to reduce the number of auxiliary variables required.  It has been shown that every self-dual code is equivalent to a graph code <cit.>. Graph codes are self-dual codes associated with a simple undirected  graph, their generator matrix can be written as $G^T = \begin{pmatrix}
I_n & B
\end{pmatrix}$, where $B$ is the adjacency matrix of the corresponding graph.

Using this representation of the generator matrix, the $i$-th term in the optimization problem (<ref>) becomes $x_i^2 + \beta_i(\beta_i -x_i)$. We can replace $\beta$  with $b$,  as defined in (<ref>), using the following identity

\begin{equation}
x_i^2 -\beta_i(\beta_i -x_i) = x_i^2 +\min_{u_i}(b_i-2u_i)(b_i-2u_i-x_i),
\end{equation}

where $u_i$ is an integer variable and can be represented in terms of binary variables as before. The corresponding $Q$ matrix then becomes:

\begin{equation}
Q = \begin{pmatrix}
I_{n}+B^2 -B  & 2\left(\boldsymbol{w_s} \otimes (I_{n}-2B)\right)^T \\
2\boldsymbol{w_s} \otimes (I_{n}-2B)  &  8\boldsymbol{w_s}\boldsymbol{w_s}^T \otimes I_n
\end{pmatrix}.
\end{equation}

§.§ Reformulating constraints

§.§.§ Trade constraints for penalty term

Different stabilizer codes are called equivalent if their stabilizer generators can be mapped to each other using conjugation by local Clifford operations[Local Clifford operations are tensor products of single-qubit unitaries
that map Pauli operators to themselves under conjugation.], and qubit permutations. The definition of code equivalence ensures that equivalent codes have the same minimum distance. Code equivalences, together with the freedom to select generators of the normalizer, can be used to bring the normalizer  matrix $G$ into various convenient forms <cit.>.

A particular convenient form to which the normalizer matrix can be efficiently  brought by a quick classical algorithm is $G = \begin{pmatrix}    L & H^T \end{pmatrix}$, where $L$ and $H^T$ are $2n \times 2k$ and $2n \times (n-k)$ binary matrices, respectively. Here, $L$ represents $2k$ linearly independent logical Pauli operators, while $H^T$ is the parity check matrix containing generators of the stabilizer group. For non-self-dual $(k \neq 0)$ codes, the constraint $(G^T\Lambda G)\boldsymbol{x} \neq \boldsymbol{0}$ can be imposed by penalizing  $\boldsymbol{x}$ with vanishing first $2k$ components $x_i$, that  correspond  to logical operations.  This is achieved by adding the following quadratic term to  the utility function

\begin{equation}
\min_m n\left(\sum_{i=1}^{2k} x_i - m\right)^2,
\quad 1\leq m \leq 2k.
\label{constraint_term}
\end{equation}

Here, $m$ is an auxiliary integer variable. Multiplying by $n$ ensures that the penalty is greater than the minimum distance (since $d < n$). To express the cost function in terms of binary variables, we trade the integer variable $m$ for $r$ auxiliary binary variables using its binary representation: $ m = 1+ \sum^{r}_{j=1} 2^{j-1} e_{j}, \ e_{j} \in \{0,1\}$. It is enough to introduce $r = \lceil \ln{2k} \rceil \leq \lceil \ln{n} \rceil $ binary variables to cover the required range for $m$.

For Self-Dual codes, the constraint $\boldsymbol{x} \neq 0^n$ can be satisfied by a term similar to (<ref>) with the sum going from $i=1$ to $n$ instead. $m$ in this case can be represented in terms of $r = \lceil \ln{n} \rceil $ binary variables. The number of auxiliary variables introduced in both cases is bounded by $\lceil \ln{n} \rceil$.

To incorporate these penalty terms into our optimization problem, we extend the  vector of variable to

\begin{equation}
\boldsymbol{z}= \left( x_1,x_2, \dots ,x_{n+k}, \ t_{1,1},t_{2,1}, \dots ,t_{n,1}, \dots ,t_{1,s},\dots, t_{n,s},u_{1,1},\dots, u_{n,s},v_{1,1},\dots, v_{n,2s}, e_1, \dots, e_r \right).
\end{equation}

The minimization problem becomes

\begin{equation}
d = n + \min_{\boldsymbol{z}} \boldsymbol{z}^T {\cal Q}\,\boldsymbol{z},
\end{equation}

with the following rows and columns added to $Q$:

\begin{equation}
{\cal Q} :=
\begin{array}{r} n+k\{ \\ \\ r \{ \end{array} \!\!\!\! \left(
\begin{array}{cc|c}
\raisebox{0ex}[1.5ex]{$\overbrace{n\left(\boldsymbol{v v^T} -2 \boldsymbol{v}\odot I_{n+k} \right)}^{n+k}$} &
0 &
\raisebox{0ex}[1.5ex]{$\overbrace{-2n \boldsymbol{v w_r^T}}^{r}$}\\
0 & 0 & 0  \\ \hline
-2n \boldsymbol{w_r v^T} & 0 & 4n(\boldsymbol{w_r w_r^T} + \boldsymbol{w_r} \odot I_r)  \end{array} \right) +
\begin{pNiceArray}[first-row, first-col]{wc{1.2cm}c|c}[margin]
&\Block{1-2}{}&  &  r \\
& \Block{2-2}{Q} & & \Block{2-1}{0} \\
&&&\\ \hline
r & \Block{1-2}{0}& & 0 \\
\end{pNiceArray}.
\end{equation}

The column vector $\boldsymbol{v}$
is defined by $v_i=1$ for $i\leq 2k$ and $v_i=0$ otherwise. By $\odot$, we denote the Hadmard (element-wise) product, $(A \odot B)_{ij} := A_{ij}B_{ij}$. Introducing the  penalty term which implements the constraint on $\boldsymbol{x}$ increases the number of binary variables only by $\ln{n}$, preserving the asymptotic  growth of the total number of variables of order $O \left(n  \ln{n/n_0}   \right)$.

§.§.§ Trade constraints for multiple cost functions

An alternative approach to incorporate the constraints is to  go through all values $i=1,\dots, 2k$, and for each value of $i$ set $x_1, \dots ,x_{i-1} = 0$ and $x_i=1$, and solve the resulting QUBO for the remaining variables. The best (minimal) value among the $2k$ QUBOs  would be the global minimum of the original constrained problem. This approach requires solving $O(n)$ distinct QUBOs of sizes $O \left(n  \ln{n/n_0}   \right)$, which is  normally less favorable than the approach presented in the previous section. However, it proves to be highly efficient for specific classes of code, particularly for circulant self-dual codes.

Circulant self-dual codes is a family of  codes defined by a normalizer matrix admitting the form $G^T = \left( I_{n}|B \right)$, with $B$ an $n \times n$ binary circulant matrix <cit.>.[Circulant matrix is a matrix with each row being a cyclic shift of the row above.] These codes are of interest since many  optimal self-dual codes discovered to date belong to this family <cit.>.[Codes that meet the linear programming upper bound on  minimum distance  are called extremal, while those that achieve the highest possible minimum distance but do not satisfy the linear programming upper bound are referred to as optimal.] To meet the constraint $\boldsymbol{x} \neq 0^{n}$, it is sufficient to set $x_1 =1$ and solve QUBO for the remaining variables.

§ ADIABATIC QUANTUM COMPUTING AND QUANTUM ANNEALING

Adiabatic Quantum Computing (AQC) is a universal model for quantum computation based on the adiabatic theorem. AQC starts by initializing the quantum processor in an easy-to-prepare ground state of some chosen initial Hamiltonian, $H(0) = H_I$. Subsequently, the system's Hamiltonian is slowly evolved over an extended time, $t_a$, transitioning towards a problem Hamiltonian, $H(t_a)=H_p$, whose ground state encodes the solution to a problem of interest. The adiabatic theorem guarantees that the evolved state $\ket{\Psi(t_a)}$ will have a large overlap with the ground state of the problem Hamiltonian provided  that the evolution was done slow enough <cit.>. In fact, this overlap, often termed the probability of success, $P_s := \abs{\braket{0,t_a}{\Psi(t_a)}}^2 $ satisfies

\begin{equation}
\lim\limits_{t_a \to \infty} P_s = 1.
\end{equation}

Here, the state $\ket{m,t}$ satisfies $H(t)\ket{m,t} = E_m(t) \ket{m,t}$ with $E_m(t)$ denoting the $m$-th lowest eigenvalue of $H(t)$ .

For Hamiltonians that can be written as $H(t)=\Tilde{H}(\gamma)$, where $\gamma:=t/t_a \in [0,1]$, and $\Tilde{H}$ has no explicit dependence on $t_a$,  the adiabatic condition can be expressed as

\begin{equation}
t_a \gg  \frac{\hbar \ \mathcal{E}}{g^2_{\text{min}}}, \quad \mathcal{E}:= \max\limits_{\gamma} \abs{\left\langle 1,\gamma \middle| \partial_{\gamma} \Tilde{H} \middle| 0,\gamma \right\rangle}, \quad g_{\text{min}} := \min\limits_{\gamma} \abs{E_1(\gamma) - E_0(\gamma)}.\end{equation}

This condition ensures that $P_s$ is finite for finite values of $t_a$ <cit.>.  The performance of AQC depends on the chosen adiabatic path, specifically the choice of the initial Hamiltonian and the trajectory of evolution towards the problem Hamiltonian. For Hamilotnians exhibiting final  ground state degeneracy, $g_{\text{min}}$ is calculated between the instantaneous ground state and the first excited state that does not evolve eventually into the final ground state subspace <cit.>. In this scenario,  we define the probability of success, $P_s$ , as the overlap between the final evolved state and the subspace of final degenerate ground-states.  Notably, AQC has been shown to be polynomially equivalent, in terms of the resource overhead, to universal gate-based quantum computation <cit.>.

Quantum Annealing (QA) is a heuristic approach that follows the same process as AQC but relaxes the condition of  adiabaticity. The departure from adiabatic evolution introduces excitations in the system's state, voiding the guarantee of convergence to the ground state of the problem Hamiltonian. To mitigate the limitations of nonadiabaticity, the algorithm repeats the entire evolution multiple times, generating various candidate final states. The algorithm then reports the best overall result. The computational power of QA and whether it offers a speedup over classical computation is currently an open area of research <cit.>, with studies comparing the performance of QA  against classical algorithms like Simulated Annealing (SA) and Simulated Quantum Annealing (SQA) <cit.>.

Quantum annealers are devices designed to run the QA algorithm. D-Wave systems have introduced multiple commercially available quantum annealers <cit.>. In this paper, we focus on the D-wave Advantage 4.1 quantum annealer, which operates by evolving a mixed-field quantum Ising Hamiltonian of the form

\begin{equation}
H(\gamma) = \underbrace{-\frac{\mathcal{A}(\gamma)}{2} \sum_i X_{(i)}}_\text{initial Hamiltonian}
+ \underbrace{\frac{\mathcal{B}(\gamma)}{2} \left( \sum_i h_i Z_{(i)} + \sum_{i>j} J_{ij} Z_{(i)} Z_{(j)}\right)}_\text{problem Hamiltonian}.
\end{equation}

Here, $\mathcal{A}(\gamma)$ and $\mathcal{B}(\gamma)$ are functions that control the transition between the initial and problem Hamiltonians, satisfying $\mathcal{A}(0)\gg \mathcal{B}(0)$ and $\mathcal{A}(1)\ll \mathcal{B}(1)$. The parameters $\gamma$ and $t_a$ are referred to as the annealing parameter and the annealing time, respectively. The existence of a map between any QUBO cost function and the problem Hamiltonian of the quantum Ising model <cit.> makes the D-Wave Advantage 4.1 quantum annealer a natural tool to solve the QUBOs constructed in the previous section.

§.§ Necessary runtime analysis

The  runtime necessary to evaluate the optimum for the worst-case instance of an NP-complete optimization problem is  believed to be exponentially large in problem size for both classical and quantum  algorithms. Hypothetically however, it could be that quantum algorithms offer a modest (sub-exponential) speedup. It is also conceivable that quantum algorithms may offer systematically better performance on average, i.e. for typical instances <cit.>. In practice there are many factors affecting   AQC runtime necessary to solve a  QUBO problem e.g. associated with the circulant self-dual optimal stabilizer codes.
We discuss the most important of them below, assuming for simplicity
that the quantum device possesses sufficient connectivity to enable direct mapping of QUBO variables onto the physical qubits. The effects of limited connectivity, which are by far the most important, will be discussed in later sections.

a) Graph of annealing time required to achieve $P_s = 0.9$, calculated using the form of $\mathcal{A}(s)$ and $\mathcal{B}(s)$ utilized by the D-Wave Advantage 4.1 quantum annealer. For larger code sizes, the dimensions of the corresponding Ising Hamiltonian exceeded  computational feasibility. b) Graph depicting the relationship between annealing time and ground-state degeneracy of the problem Hamiltonian.

It is well-known that the degeneracy of the ground state, even in absence of quantum noise and decoherence, may significantly influence necessary annealing time $t_a$.
Given that circulant codes are cyclic, our problem Hamiltonian  $H(\gamma=1)$ always exhibits multiple ground states. Ground state degeneracy can both increase and decrease $t_a$, depending on the degeneracy of the excited states <cit.>. To address this question in our case we numerically simulated the time-dependent Schrodinger equation describing the annealing process for codes of sizes $n=3-8$. In these simulations annealing time was varied until the success rate $P_s= 0.9$ was achieved. The results, shown in the panel a) of Fig. 1 are consistent with the inverse correlation between the ground state degeneracy and minimal necessary $t_a$.

This conclusion is very preliminary, as the degeneracy of final Hamiltonian is merely one of the factors. Another important factor is the smallest value of gap of $H(\gamma)$, which may change significantly for large $n$ (we checked that for codes of sizes $n=3-7$ it does not). For the Hamiltonian employed, it can be demonstrated that $\mathcal{E}$ scales polynomially with the problem size, as follows from the observation that the eigenvalues of $H_I$ and $H_P$ are also at most polynomial. Yet,  understanding the scaling of $g_{\text{mim}}$ with $n$ is usually a very difficult problem, see <cit.>.

Our simulations shown in Fig. 1, panel b) suggest that theoretically required  annealing times for the small code sizes $n\sim 3-8$  fall within the nanosecond range. In practice, as we discussed in the next section, reasonable necessary annealing time on the D-Wave Advantage 4.1 platform for these problems is of the order of microseconds. This suggests
theoretically quantum annealer is operating firmly within the adiabatic regime, but in practice its quality of performance is dominated by other factors related to quantum noise. The factors discussed above are much less important, at least for moderate problem sizes $n$. This instructs our approach in the next section to optimize the annealing time $t_a$ empirically.

§ RESULTS OF SIMULATIONS ON THE D-WAVE QA PLATFORM

§.§ Quantum-only implementation

We tested our QUBO-based algorithm by running implementations on the D-Wave Advantage 4.1 quantum annealer. We ran instances of optimal circulant self-dual codes of small length, taken from the database composed by M. Grassl <cit.>. The constraint $\boldsymbol{x} \neq 0^n$  was imposed by setting $x_1 =1$ as described in section <ref>. To quantify how close the quantum annealer got to the optimal value, we calculated the so-called approximation ratio for various code sizes. The approximation ratio is the ratio between the smallest value generated by the quantum annealer and the actual minimum distance. As shown in Figure <ref>, the quantum annealer was able to successfully find exact minimum distance for codes of sizes up to $12$. For larger code sizes, despite numerous re-runs, QA could only find conservative upper bounds. Our results strongly suggest that QA alone can not currently compete with purely classical algorithms/implementations, but there could be a room to use QA  as a tool to obtain rigorous upper bounds. Efficiency of this method should be addressed separately by comparing the upper bounds with the theoretically known ones <cit.>.

Approximation Ratio (AR) of best output generated by runs of quantum annealer solving the QUBO associated with optimal circulant codes of length $6\leq n\leq 15$. AR equals one means the algorithm found the global minimum (exact value of the minimum code distance). Larger AR means that despite numerous re-runs, the algorithm failed to find the globabl minimum and only produced a conservative bound.  We performed 40 runs (one run is a 100 anneals) for each code size.

To systematically analyze the performance of quantum annealer,
we study how the quality of outcome depends on the annealing time $t_a$,
a free parameter of implementation of the quantum algorithm. At least naively, longer annealing time would make time evolution more adiabatic. But since the platform is noisy, a longer annealing time can lead to an accumulation of errors. Thus, we can expect that certain values of $t_a$ would yield optimal performance.

To study $t_a$ dependence, we plot the rate of success for the QUBO associated with the given code of length $n=9$ using various annealing times. The rate of success is defined as the proportion of runs where the optimal result was obtained out of the total number of runs conducted. As shown in Figure <ref>, the rate of success initially increases as we increase the annealing time, reaching its peak at approximately $4.5 \ \mu s$. However, increasing the annealing time beyond this point resulted in deteriorating performance. Considering the results from the previous section, which suggested that, at least theoretically, we are operating within the adiabatic regime, this deterioration in performance is likely due to accumulated noise.  Empirically, we find that an annealing time of $t_a = 4.5 \ \mu s$ is close to optimum for most code sizes, and therefore, we  use this value in all of our simulations.

The success rate of the QUBO-based algorithm run on D-Wave Advantage 4.1 as a function of different annealing times $t_a$. Each data point represents the success rate, calculated as the number of times the optimal value was achieved out of 40 runs (each comprising a 100 anneals). All runs were optimizing the same QUBO associated with an optimal circulant code of length $n=9$. The annealing times ranged from $0.5$ to $8.5  \mu s$ in increments of $1 \mu s$.

Suboptimal performance of quantum annealer for relatively small code sizes $n\geq 13$ can likely be attributed to  limited connectivity between physical qubits of the platform. If two physical qubits are connected, the platform can directly implement two-qubit gates.
Thus, for best performance,  any two binary variables in the cost function that share a cross-term should be mapped to two physically connected qubits. However, in many practical scenarios, each QUBO variable can share cross-terms with numerous other variables, exceeding the connectivity limitations of individual device qubits in the contemporary available platforms. A common workaround is to map a single QUBO variable to a chain of physical qubits, a process referred to as embedding <cit.>.  This leads to an overhead in the number of qubits employed by the platform. Figure <ref> shows that the number of necessary physical qubits increases approximately quadratically with the size of the original QUBO problem, associated with the code. The increase in the number of necessary qubits significantly complicates successful implementation of the QA algorithm. As the number of physical qubits grows, the corresponding search space also expands. Furthermore, since current quantum devices are noisy, each additional qubit introduces additional sources of noise that can affect the final result of the algorithm.

Mean number of physical qubits employed by the QA platform (i.e. number of qubits after the so-called embedding) vs. the number of binary variables of the underlying QUBO problem.
Standard deviation, due to heaurstic nature of the embedding algorithm, is shown as the error bars. It is usually small.
To generate the plot, we used
optimal circulant codes of different sizes, reformulated the minumum distance problem as QUBO and then embed it onto the QA platform 10 times. We note, that QBUO associated with different codes can have different  connectivity requirements, which can affect the number of qubits after embedding.
The data points (blue) conform well to quadratic fit (orange), which suggests that on average the number of qubits after embedding grows quadratically.

§.§ Hybrid implementation

To mitigate challenges of limited connectivity, we shift our focus to hybrid quantum-classical implementation. Specifically, we use Qbsolv algorithm, a hybrid algorithm developed by the D-Wave Systems. Qbsolv employs a classical algorithm to partition QUBO problems into smaller and more manageable pieces. These “subproblems”  can be solved using QA, or any other suitable QUBO solving algorithm, before stitching them together <cit.>. This approach offers two significant benefits. First, it allows solving problems of bigger size, exceeding current limits of  quantum platforms. Second, the subproblems have fewer connectivity constraints, which reduces the number of necessary physical qubits,  and in this way decreases the source of errors.

To compare performance of hybrid approach against a fully classical one, we conducted runs solving subproblems using purely classical Simulated Annealing (SA) algorithm as well as QA using the D-Wave platform. In what follows we refer to these runs as  “SA-assisted” and “QA-assisted.” To ensure a fair comparison,  for each run we imposed  the same overall time limit, which reflects the time-cost from the end-user point of view.  However, it should be noted that the QA-assisted runs experience hidden latency of unknown duration, due to remote execution of QA on the D-Wave platform. In contrast, SA-assisted runs were executed locally and did not experience this issue. Therefore, although end-user time was the same, the QA implementation was effectively subject to a more stringent time constraint.

We performed a different number of runs for different codes, but each QUBO problem, associated with a given code, was minimized the same number of times using QA and then SA-assisted implementations. Since the subproblems are of the similar size to the problems discussed in the previous section, we use an annealing time of $t_a = 4.5 \ \mu s$ also for the hybrid quantum-assisted simulations.

Our results, shown in Figure <ref>, demonstrate a comparable performance between the two approaches. For codes of length up to $n=20$, both approaches managed to find exact minimum distance $d$. For longer codes of length $n\lesssim 30$, both algorithms, with some exceptions, typically find the exact $d$, but as $n$ increased, there were codes for which either implementation failed despite many repetitive runs. Of the $17$ codes for which either of the implementations did not find the exact result, in approximately $58\%$ of cases both algorithms yielded only conservative bounds. This may suggest that for certain codes, the resulting QUBO problems are inherently hard, regardless of which implementation is used.
In certain cases, for codes of length $n=21, 26, 40$ and $49$ ($24\%$ of all cases when either approach failed), the hybrid approach succeeded in finding the QUBO global minimum while the classical approach failed. On the contrary, the reverse scenario occurred only three times ($18\%$ of all cases), for codes of length $n=29,31$ and $38$.

Approximation Ratio (AR) for minimum distance $d$ for codes of varying length, obtained with help of QUBO formulation, optimized using hybrid QA-assisted or purely classical SA-assisted approach. For each code length $n$, the number of QA and SA-assisted runs were the same (but different for different $n$).
Due to limited computational resources, We bypassed certain code lengths in favor of larger values.

§ DISCUSSION

In this work we introduced an approach to reformulate the problem of finding the minimum distance of a binary quantum stabilizer code as a QUBO (quadratic unconstrained binary optimization) problem. The underlying motivation behind this was to ultimately harness the capabilities of contemporary quantum platforms to help design better quantum error-correcting protocols. As a result, we reformulated a well-known NP-hard problem of coding theory in a way suitable for solving using quantum processors. We then investigated the practical advantages of this approach by running problem instances on the D-Wave Advantage 4.1 platform. Our QUBO reformulation requires  $O(n\ln\, n)$ optimization variables, representing only a logarithmic multiplicative overhead.

The QUBO formulation of the minimum distance problem offers a conceptual and possibly practical advantage. First, QUBO is an active area of research. Therefore practical implementation of our QUBO-based approach can benefit from many existing heuristics and speedups. Second, the QUBO-based formulation can be solved using quantum, classical, or hybrid methods. While our study was focused on investigating practical advantage of solving QUBO using quantum annealing, one can envision in the future using coherent gate-based quantum devices and  algorithms like the Quantum Approximate Optimization Algorithm (QAOA) <cit.> to tackle the minimum distance and similar problems.

In our simulations, we compared purely quantum, classical, and hybrid approaches to solving the QUBO associated with the minimal distance problem. Our focus was on optimal self-dual codes of length $n$,
i.e. self-dual codes with the largest possible minimal distance $d$ for given $n$. Minimum distances for theses codes with small or moderate $n$ are known, providing a good testing ground. We found that the purely quantum approach fails to produce the optimum (only providing conservative values for $d$), already for $n\gtrsim 10$. Purely classical and quantum-classical hybrid approaches have shown much better and mutually comparable performance, yielding correct values of $d$ for codes up to $n\sim 60$, see Figure <ref>. In the cases where the algorithms failed to produce the correct distance, they provided rigorous upper bounds. This suggests that for even larger $n$, when deterministic approaches become completely impractical, heuristics can be employed to solve these QUBOs to find rigorous upper bounds on the minimum distance of specific code instances. Whether  these upper bounds can be made tighter than theoretical upper bounds established by linear programming techniques <cit.> remains an open question.

Overall, our simulations have demonstrated potential of hybrid approach, that combines quantum annealing with classical algorithm used to partition the problem into smaller pieces. Our study highlights that the limited connectivity of quantum devices, along with the associated embedding problem, represents a primary obstacle on the path to practical advantage of QA.  As is shown in Figure <ref>, embedding the QUBOs onto the quantum platform  introduces roughly a quadratic overhead in the number of variables. Due to intrinsic noise, the overhead in the number of qubits not only enlarges the search space but also leads to significantly more errors, which results in the suboptimal minimum.

Running simulations on the D-Wave Advantage 4.1 platform, we concluded that the QUBO-based approach currently does not offer a practical advantage over the classical deterministic minimum distance  algorithm introduced by G. White and M. Grassl <cit.>. However, as the  fault-tolerant quantum devices continue to develop, this situation might change. While classical deterministic algorithm would require the runtime exponentially growing with $n$, AQC-based approach may
show better scaling. Understanding this scaling theoretically, and thus certifying possible quantum advantage of the QUBO-based approach, would be an important task for the future.

To conclude, we proposed a reformulation of a theoretically and practically important NP-hard problem as a QUBO problem, contributing to an ongoing effort to reformulate various problems in a form suitable for quantum platforms. Looking ahead, an important step would be to extend our approach to other problems of classical and quantum information theory and discrete mathematics. In particular, the shortest lattice vector problem, a known NP-hard problem, which is conceptually similar to code minimum distance problem discussed in this work, would be a natural  next candidate.

§ ACKNOWLEDGMENTS

We thank D-Wave for providing access to their quantum cloud services, Leap. This research is supported by the NSF grant PHY 2310426.