# Breakdown of the Hund’s Rule in CuFeAs

Ze-Yi Song  Shanghai Key Laboratory of Special Artificial Microstructure
Materials and Technology, School of Physics Science and engineering, Tongji
University, Shanghai 200092, P. R. China  Xiu-Cai Jiang  Shanghai Key
Laboratory of Special Artificial Microstructure Materials and Technology,
School of Physics Science and engineering, Tongji University, Shanghai 200092,
P. R. China  Xiao-Fang Ouyang  School of Electronic and Electrical
Engineering, Shangqiu Normal University, Shangqiu 476000, P. R. China  Yu-
Zhong Zhang  Email<EMAIL_ADDRESS>Shanghai Key Laboratory of Special
Artificial Microstructure Materials and Technology, School of Physics Science
and engineering, Tongji University, Shanghai 200092, P. R. China

###### Abstract

The ground-state properties of CuFeAs were investigated by applying density
functional theory calculations within generalized gradient approximation (GGA)
and GGA+U. We find that the bicollinear antiferromagnetic state with
antiparallel orbital magnetic moments on each iron which violates the Hund’s
rule is favored by the on-site Coulomb interaction, which is further
stabilized by Cu vacancy. The magnetic ground state can be used to understand
weak antiferromagnetism in CuFeAs observed experimentally. We argue that
breakdown of the Hund’s rule may be the possible origin for reduced magnetism
in iron pnictides, rather than magnetic fluctuations induced by electronic
correlations.

##  I introduction

The discovery of the high-temperature superconductivity in iron-based
compounds with transition temperature $T_{C}$ up to 55 K have attracted
tremendous research interest for the last decade Kamihara  _et al._ (2008);
Rotter  _et al._ (2008); Zhi-An  _et al._ (2008); Stewart (2011). Similar to
the phase diagram in cuprate Keimer  _et al._ (2015), the superconductivity in
iron-based superconductors always emerges in close proximity to a state with
antiferromagnetism de la Cruz  _et al._ (2008); Basov and Chubukov (2011);
Scalapino (2012); de’ Medici  _et al._ (2014). Since previous theoretical
studies have ruled out the phonon-mediated pairing Boeri  _et al._ (2008), it
is widely believed that the superconductivity in iron-based superconductors is
unconventional and has a magnetic origin.

Therefore, many efforts have been spent on revealing the origin of magnetism
in iron-based superconductors. As continuous debates between itinerant
scenario Mazin  _et al._ (2008); Dong  _et al._ (2008); Mazin (2010) and
localized picture Yildirim (2008); Si and Abrahams (2008) remain unsettled, a
compromising explanation arised where coupling of itinerant electrons and
localized spins is taken into account Johannes and Mazin (2009); You  _et al._
(2011). In spite of the above open discussions, two consensuses have been
reached. One is that density functional theory (DFT) calculations can
qualitatively describe the magnetic properties of parent states in both iron
pnictides and iron chalcogenides, although the magnetic moments are always
overestimated due to the fact that the intermediate strength of electronic
correlations in open $d$ shell of iron atom can not be properly captured in
the approximation of the functionals. The other is that, in contrast to the
cuprates, the effect of electronic correlations is adopted to come from the
Hund’s rule coupling Haule and Kotliar (2009); Yin  _et al._ (2011); Lanatà
_et al._ (2013), rather than on-site Coulomb repulsion. And the metallic
states of iron-based superconductors are called Hund’s metal where the Hund’s
rule of maximum multiplicity is supposed to be valid Georges  _et al._ (2013).

However, above consensuses have been seriously challenged since the discovery
of CuFeAs. The magnetic susceptibility measurements showed that it is
antiferromagnetic with Néel temperature $T_{N}$ of around 9 K Thakur  _et al._
(2014); Li  _et al._ (2018). The antiferromagnetism was further demonstrated
by neutron diffraction experiments Zou  _et al._ (2017); Kamusella  _et al._
(2017), where either an unusual G-type antiferromagnetic order or proximity to
an antiferromagnetic instability was proposed. Though ferromagnetism was also
reported in the literature Qian  _et al._ (2015); Lv (2009), it was pointed
out that the weak ferromagnetism probably comes from a ferromagnetic component
of a canted antiferromagnetic state Zou  _et al._ (2017); Qian  _et al._
(2015). Therefore, while the type of antiferromagnetic order is still unclear,
it seems conclusive that the ground state of CuFeAs is antiferromagnetic
experimentally. Nonetheless, an early theoretical study based on DFT
calculations suggested that this compound is a ferromagnet Wang  _et al._
(2016), in stark contrast to the experimental observations Thakur  _et al._
(2014); Li  _et al._ (2018); Zou  _et al._ (2017); Kamusella  _et al._ (2017),
which casts doubts on the existing consensus. Obviously, further studies based
on DFT calculations are required to resolve the contradiction, and to clarify
the magnetic structure in CuFeAs, as well as to verify the role of Hund’s
coupling.

Here, the nature of magnetism of CuFeAs is investigated by applying DFT
calculations. We find that the ground-state magnetic structure of CuFeAs are
controlled by As height $h_{\text{As}}$ from iron plane, similar to other
iron-based superconductor parent compounds, and the critical height
$h_{\text{c}}$ of 1.612 Å is identified. If $h_{\text{As}}<h_{\text{c}}$, the
ground state is in a collinear antiferromagnetic (CAFM) state. On the
contrary, when As height is larger than $h_{c}$, the on-site Coulomb
interaction is crucial to involve in order to correctly account for observed
antiferromagnetic state. It is found that bicollinear antiferromagnetic (BAFM)
order gives lowest total energy among the states we studied, which becomes
even more favorable in the intermediate value of on-site Coulomb interaction
after introducing Cu vacancy and shows weak ferrimagnetism where total
magnetic moment turns out to be nonzero. The small magnetic moment per iron is
ascribed to the violation of Hund’s rule
[Duetothepresenceoftetragonalcrystalfieldiniron-basedcompounds\text{;
}togetherwiththenegligiblestrengthofspin-orbitcouplinginiron\text{;
}theorbitalangularmomentumof$3d$orbitalsiscompletelyquenchedinCuFeAs.Therefore\text{;
}thesecondandthethirdHund'srulesareirrelevantinCuFeAs\text{;
}andthroughoutthepaper\text{; }theHund'srulereferstothefirstone]Hundsrule
where antiparallel orbital magnetic moments on each iron are present. Our
results can be applied to fully understand the experimental results Thakur
_et al._ (2014); Li  _et al._ (2018); Zou  _et al._ (2017); Kamusella  _et
al._ (2017); Qian  _et al._ (2015).

CuFeAs is isostructural to 111-type iron-pnictide superconductor parent
compounds LiFeAs Tapp  _et al._ (2008) and NaFeAs Li  _et al._ (2009a). It is
characterized by large $h_{\text{As}}$ in comparison to other iron pnictides,
varying from $h_{\text{Kam}}$=1.53 Å Kamusella  _et al._ (2017) to
$h_{\text{Li}}$=1.57 Å Li  _et al._ (2018) and $h_{\text{Thakur}}$=1.74 Å
Thakur  _et al._ (2014), finally to $h_{\text{Zou}}$=1.80 Å Qian  _et al._
(2015). Moreover, it was reported to be nonstoichiometric Thakur  _et al._
(2014); Qian  _et al._ (2015); Zou  _et al._ (2017); Kamusella  _et al._
(2017); Li  _et al._ (2018), namely the Cu vacancies are always present.
Finally, like other iron-based superconductor compounds, CuFeAs is a material
with large Sommerfeld coefficient, indicating that electronic interactions are
not negligible Qian  _et al._ (2015).

##  II method

The DFT calculations were performed using full-potential linearized augmented
plane wave method as implemented in the Wien2k code Schwarz  _et al._ (2002).
We adopted generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhhor
Perdew  _et al._ (1996) for the exchange-correlation potentials. In order to
determine the magnetic ground state, we have studied nonmagnetic (NM),
ferromagnetic (FM), and three distinct antiferromagnetic configurations,
including Néel antiferromagnetic (NAFM), CAFM, and BAFM orders Yin  _et al._
(2010); [SeealsocartoonsforthesemagneticstatesintheappendixA]cartoon. We
employed $\sqrt{2}\times\sqrt{2}\times 1$ and $2\times 1\times{1}$ unit cells
for CAFM and BAFM states, and primitive cell for the other states,
respectively. The Brillouin zone integration is carried out with a k mesh of
$24\times{24}\times{18}$ for NM, FM, and NAFM states, and
$24\times{24}\times{20}$ for CAFM phase as well as $16\times 32\times 20$ for
BAFM phase, respectively. Furthermore, the influence of on-site Coulomb
interaction on magnetic stabilities of CuFeAs was investigated by GGA+U
approach Liechtenstein  _et al._ (1995). The around mean-field double counting
Czyżyk and Sawatzky (1994) is employed as CuFeAs is a correlated metal. Unless
specified otherwise, the Hund’s coupling $J=U/10$ Miyake  _et al._ (2010) was
used throughout all our GGA+U calculations [Theexplicitinclusionofon-
siteCoulombrepulsion$U$andexchangeparameter$J$whichfavorsthefirstHund'sruleandthusisso-
calledtheHund'scouplingtermallowsustoinvestigatemulti-
orbitalphysicsinCuFeAs.IncontrasttotheconventionalunderstandingwheremaximummultiplicityshouldappearinthepresenceoftheHund'scoupling\text{;
}breakdownoftheHund'sruleoccursinCuFeAs\text{;
}resultinginanexoticBAFMstatewithsmallmagneticmoment]methods. The conclusions
remain valid if $J$ is fixed at $U/4$. Here, the experimental lattice
constants reported by Thakur et al. Thakur  _et al._ (2014) was used and also,
the conclusions will not be altered when other experimental lattice constants
are used. The $x$ ($y$) of local coordinate system for each iron is along the
closest Fe-Fe bond direction.

Figure 1: (Color online) (a) The calculated total energies of bicollinear
antiferromagnetic (BAFM), collinear antiferromagnetic (CAFM), Néel
antiferromagnetic (NAFM), ferromagnetic (FM) and nomagnetic (NM) states as a
function of As height $h_{\text{As}}$. Here, the energies of NM state were set
to zero. (b) By mapping the energies of ferromagnetic and various
antiferromagnetic states onto a Heisenberg model, the nearest-, next-nearest-,
and third-nearest-neighbor exchange couplings $J_{1}$, $J_{2}$, $J_{3}$ were
derived. The inset denotes the cartoon of $J_{1}$, $J_{2}$, $J_{3}$ exchange
interactions.

##  III results

Since it is well known that the magnetic properties of parent states of iron-
based superconductors are sensitive to As height $h_{\text{As}}$ measured from
iron plane Yildirim (2009); Yin  _et al._ (2008), we first investigated the
effect of $h_{\text{As}}$ on ground-state magnetic structure by only varying
the As $z$ positions and leaving all the other internal coordinates unchanged.
Fig. 1 (a) shows the calculated total energies of different magnetic
configurations as a function of As height within GGA. The energies of the NM
state were set to zero. Starting from As height $h_{\text{Kam}}$ of 1.53 Å,
CAFM state gives lowest total energy. When $h_{\text{As}}$ exceeds critical
height $h_{c}$ of 1.612 Å, FM state, rather than any AFM states, becomes
stablest. The appearance of phase transition from a CAFM state to a FM one
suggests that the magnetism is strongly dependent on $h_{\text{As}}$ in
CuFeAs, similar to the other parent compounds of iron-based superconductors
Yin  _et al._ (2008); Yildirim (2009); Moon and Choi (2010).

However, remarkable conflict can be found between the calculated results and
the experimental data. As can been seen in Fig. 1 (a) where various As heights
reported in different experiments are shown, at $h_{\text{As}}$ less than
$h_{\text{c}}$, our theoretical calculations point to a CAFM ground state,
which is mainly consistent with the experimental results where either long-
range Li  _et al._ (2018) or short-range Kamusella  _et al._ (2017) AF order
was observed, though the magnetic structure has not been determined
experimentally yet. But, in the cases of $h_{\text{As}}$ greater than $h_{c}$,
while the theoretical ground state is strongly FM, it was inferred from
experiments Thakur  _et al._ (2014); Qian  _et al._ (2015); Zou  _et al._
(2017) that CuFeAs should be an antiferromagnet probably with ferromagnetic
components. The contradiction raises a great challenge to the existing theory
Moon and Choi (2010); Yin  _et al._ (2010), which seems valid for the whole
family of iron-based superconductors including both iron pnictides and iron
chalcogenides, where various magnetic ground states can be correctly accounted
for after different anion height is considered. In fact, the theory of anion
height is also true for a sister compound of CuFeAs, so called CuFeSb, where
DFT calculations Wang  _et al._ (2016) and experiments Qian  _et al._ (2012);
Sirohi  _et al._ (2016) both obtained a FM ground state, which can be
attributed to a very high anion height of $>1.84$ Å. Therefore, it is urgent
to understand why CuFeAs with intermediate value of anion height is so
extraordinary that DFT calculations can not agree with experimental findings
even qualitatively.

Till now the magnetism of iron pnictides can be explained by both the local
moment picture Yildirim (2008); Si and Abrahams (2008) and the weak-coupling
scenario Mazin  _et al._ (2008); Dong  _et al._ (2008); Mazin (2010). From
local moment picture, the magnetic ground state can be effectively described
by the frustrated Heisenberg model $H=\sum_{ij}J_{ij}S_{i}S_{j}$ where
$J_{ij}$ is the superexchange interactions between local Fe moments with spin
$S_{i}$. From above data based on the DFT calculations within GGA, the
nearest, next-nearest, and next-next-nearest neighbor exchange couplings
$J_{1}$, $J_{2}$, and $J_{3}$ can be derived from the energy differences among
various magnetic states and are summarized in Fig. 1 (b). As can be seen,
$J_{1}$ remains FM, thus favors a FM order. It is strongly enhanced as a
function of increased $h_{\text{As}}$. While $J_{2}$ is AFM and plays a
dominant role at small $h_{\text{As}}$, it drastically reduced in the vicinity
of $h_{c}$ and eventually turns into FM nature at larger $h_{\text{As}}$.
$J_{3}$ is also AFM and barely dependent on $h_{\text{As}}$. If only FM, NAFM,
CAFM, and BAFM states are taken into account in the classic limit of
Heisenberg $J_{1}-J_{2}-J_{3}$ model, CAFM is energetically favorable over
other magnetic configurations when $J_{2}>-J_{1}/2$ and $J_{2}>2J_{3}$. The
conditions are satisfied when $h_{\text{As}}$ is smaller than $h_{\text{c}}$,
resulting in an agreement between our first-principles results and recent
experimental observations Kamusella  _et al._ (2017); Li  _et al._ (2018).
However, when $h_{\text{As}}$ is greater than $h_{\text{c}}$, the FM state
becomes the ground state owing to $J_{1}<-2J_{2}$ and $J_{1}<-J_{2}-2J_{3}$.
Since antiferromagnetism was observed in experiments Thakur  _et al._ (2014);
Qian  _et al._ (2015); Zou  _et al._ (2017) when
$h_{\text{c}}<h_{\text{As}}<1.84$ Å, it indicates that local moment picture
fails to account for the magnetic ground state of CuFeAs in the intermediate
region of anion height.

Figure 2: (Color online) the orbitally-resolved Pauli susceptibilities of Fe
3$d$ orbitals for As height equal to 1.74 Å (a) and 1.53 Å (b) along high
symmetry path of $\Gamma-\text{X}-\text{M}-\Gamma$.

In order to know if the magnetism can be understood from the weak-coupling
limit where the Fermi surface nesting plays roles, we have calculated the
orbitally resolved Pauli susceptibility $\chi^{\tau\tau}_{\tau\tau}(q)$, as
defined in Ref. Graser  _et al._ (2009); Ding  _et al._ (2013), to quantify
the nesting property. Owing to the fact that the magnetism is mainly
controlled by intraorbital particle-hole excitations, we only show in Fig. 2
the intraorbital components of the Pauli susceptibility along the path of
$\Gamma-\text{X}-\text{M}-\Gamma$ in the Brillouin zone. It is found that the
susceptibility of $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals are much smaller
than those of $d_{xy}$ and $d_{yz/xz}$ orbitals, suggesting that the magnetic
instabilities in CuFeAs are mainly contributed from the latter three orbitals.
As is well known, a prominent peak present in the Pauli susceptibility denotes
a tendency towards a certain long-range magnetic ordered state whose magnetic
configuration is determined by the position of the peak in the momentum space.
However, in the intermediate region of anion height, for instance,
$h_{\text{Thakur}}$=1.74 Å as shown in Fig. 2 (a), neither the
susceptibilities of $d_{xy}$ orbital nor those of $d_{xz/yz}$ orbitals show
any pronounced peaks. On the contrary, the plateau appearing in the
susceptibilities of $d_{xy}$ orbital along the $\text{X}-\text{M}$ path may
indicate that CuFeAs is highly magnetically frustrated due to the competitions
among enormous instabilities.

For comparison, we have calculated the Pauli susceptibility at low As height
of 1.53 Å, as shown in Fig. 2 (b). In contrast to the featureless orbital-
resolved susceptibilities at $h_{\text{As}}=$ 1.74 Å, the counterpart at low
As height shows pronounced instabilities in $d_{xy}$ orbital around wave
vector ($\pi,\pi$), indicating a strong tendency towards a antiferromagnetic
state. This is consistent with our total-energies calculations and previous
experimental observations Kamusella  _et al._ (2017); Li  _et al._ (2018). It
suggests that the magnetism in CuFeAs for $h_{\text{As}}<h_{c}$ case can be
explained by both the Fermi surface nesting scenario Mazin  _et al._ (2008);
Dong  _et al._ (2008); Mazin (2010) and the local moment picture Yildirim
(2008); Si and Abrahams (2008), similar to all the other iron pnictides.
However, if $h_{c}<h_{\text{As}}<1.84$ Å, neither the local moment picture nor
the Fermi surface nesting scenario could be applied to understand the
magnetism in CuFeAs.

Figure 3: (Color online) The effect of on-site Coulomb interaction $U$ on the
magnetic ordering in CuFeAs (a) and Cu0.5FeAs (b). As $U$ increases, the phase
transiton from FM state to BAFM state occurs at around 4.1 eV for CuFeAs (a)
and at about 3.5 eV for Cu0.5FeAs (b), respectively.

Considering the orbital degrees of freedom existing in CuFeAs, inclusion of
on-site Coulomb interactions may strongly change the spin and charge
populations among different orbitals. Therefore, in the following, we will
apply GGA+U Czyżyk and Sawatzky (1994) to allow for multi-orbital effects and
unravel the truth for magnetic ground state in CuFeAs. Fig. 3 (a) shows the
evolution of total energies of different magnetic configurations as a function
of $U$. At small $U$, the FM state gives the lowest total energy. As $U$
becomes large, the BAFM state becomes stabler than other magnetic states. The
FM-BAFM phase transition takes place at critical point $U_{c}$ of around $4.1$
eV, which is slightly smaller than the on-site Coulomb interaction $U\approx
4.5$ eV estimated from constrained local density approximation which is
comparable to those in other iron pnictides Miyake  _et al._ (2010). This
implies that the material is in the BAFM state and close to the phase boundary
between FM and BAFM states.

Furthermore, as was observed in experiments Qian  _et al._ (2015); Zou  _et
al._ (2017), CuFeAs is nonstoichiometric with Cu sites being partially vacant.
Cu vacancy will cause heavy hole doping and may alter crystal field splitting
due to the deficiency of cations in certain positions. Here, the effect of Cu
vacancy on the magnetic properties of ground state was considered by simply
removing one Cu from the primitive cell, leading to chemical formula of
Cu0.5FeAs. The total energies of various magnetic configurations as a function
of $U$ are displayed in Fig. 3 (b). Similar results are obtained as those for
the stoichiometric case, except that the critical value of $U$ is considerably
suppressed from around $4.1$ eV to about $3.5$ eV, suggesting that the BAFM
state is further stabilized when Cu vacancies are introduced.

Figure 4: (color online) (a) The averaged magnetic moment on Fe atoms for
CuFeAs and Cu0.5FeAs. The inset is the total Fe magnetic moments for Cu0.5FeAs
in the BAFM state. The orbitally-resolved magnetic moments of Fe 3$d$ orbitals
for CuFeAs (b) and Cu0.5FeAs (c). For Cu0.5FeAs in the BAFM state (c), the
averaged orbital magnetic moments over irons of two sublattices are different,
leading to the formation of ferrimagnet.

To gain deep insight into the effect of on-site Coulomb interaction and Cu
vacancies on the ground state, we have calculated total and averaged magnetic
moments on iron as well as orbitally resolved magnetic moments of 3$d$
orbitals on iron for both CuFeAs and Cu0.5FeAs. Fig. 4 (a) displays the
averaged Fe moments. It is found that the averaged Fe moment decreases as $U$
increases and at critical point $U_{c}$, a sharp drop appears. The presence of
Cu vacancy gives rise to a decrease of the averaged Fe moments in the FM state
while an increase in the BAFM state, indicating that the Cu vacancy
energetically stabilizes the BAFM state and enhances the spin polarizations in
the BAFM state. In addition, if Cu vacancy is present, the total magnetic
moments of Fe are finite in the BAFM state, as depicted in the inset of Fig. 4
(a) which is consistent with the nonzero spontaneous magnetic moment observed
by experiments Thakur  _et al._ (2014); Qian  _et al._ (2015).

Fig. 4 (b) and (c) shows the averaged magnetic moments of Fe 3$d$ orbitals for
CuFeAs and Cu0.5FeAs, respectively. The orbitally-resolved magnetic moments
are obtained by constructing atomic projectors as implemented in Wien2k code
Schwarz  _et al._ (2002). It is found that, while the magnetic moments of
$d_{xy}$ and $d_{yz/xz}$ orbitals which dominate the states close to the Fermi
level are changed slightly in FM state, those of $d_{z^{2}}$ and
$d_{x^{2}-y^{2}}$ decrease remarkably and finally become antiparallel to the
magnetic moments of $d_{xy}$ and $d_{yz/xz}$ orbitals. Such kind of breakdown
of the Hund’s rule was proposed to be a possible origin for the low
magnetization in LaFeAsO Cricchio  _et al._ (2010); Bascones  _et al._ (2010);
Liu  _et al._ (2012), where DFT calculations Cricchio  _et al._ (2010)
suggested that opposite magnetization among different orbitals is stabilized
against the Hund’s rule by the formation of large multipoles of the spin
density, while model calculations based on a five-band Hubbard model Bascones
_et al._ (2010) and a two-orbital Heisenberg model Liu  _et al._ (2012) found
that interorbital exchange interaction overrides the Hund’s rule. From our DFT
calculations and corresponding derivations of tight-binding model parameters
by Wannier90 Mostofi  _et al._ (2008), we conclude that both magnetic
multipoles Cricchio  _et al._ (2010) and interorbital hoppings Bascones  _et
al._ (2010) which induce interorbital exchanges as the Hubbard interaction $U$
is involved play roles in the formation of antiparallel magnetic moments among
five 3d orbitals on each iron in CuFeAs. Further increasing $U$, a phase
transition from the FM state to a BAFM state occurs with a significant
reduction of magnetic moment in each orbital. However, the breakdown of the
Hund’s rule remains. The weak antiferromagnetism agrees well with experimental
results Thakur  _et al._ (2014); Zou  _et al._ (2017); Kamusella  _et al._
(2017); Li  _et al._ (2018). Note in the presence of Cu vacancy, the orbital
magnetic moments on different magnetic sublattices are strongly deviated from
each other (Fig. 4 (c)), leading to a finite total magnetic moment as shown in
the inset of Fig. 4 (a).

##  IV discussions

From above investigations, it was shown that the magnetic properties of CuFeAs
vary from a CAFM state to a BAFM phase as a function of the experimental As
height measured from iron plane, where the critical height for the CAFM-BAFM
transition is $h_{\mathrm{c}}\approx 1.612$ Å. At the As height of
$h_{\mathrm{As}}>1.84$ Å, the pnictogen-height driven FM phase is expected to
be stablest Qian  _et al._ (2012); Sirohi  _et al._ (2016). This can account
for various experimental observations in CuFeAs, where all magnetic
susceptibility measurements at low external magnetic fields Li  _et al._
(2018); Thakur  _et al._ (2014); Qian  _et al._ (2015) suggest this material
is in an antiferromagnetic state for $1.57$ Å $\leq{}h_{\mathrm{As}}\leq 1.80$
Å. At $h_{\mathrm{As}}=1.53$ Å, the short-range CAFM order was observed by the
M$\ddot{\rm{o}}$ssbauer spectroscpy and muon spin resonance experiments
Kamusella  _et al._ (2017), which is also qualitatively consistent with our
results (shown in Fig. 1 (a)) since the quantum fluctuations are completely
frozen in the DFT calculations. Aware of a fact that the interlayer spin
exchange hardly affects the magnetic structure in the Fe2As2 plane for CuFeAs,
the G-type antiferromagnetic state, a three dimensional NAFM state proposed in
recent study Zou  _et al._ (2017), can not become the ground-state magnetic
ordering because the energy of NAFM state is higher than those of the other
antiferromagnetic states we considered, as displayed in Fig. 1 (a) and Fig. 3.
Besides, the agreement between our results and experiments indicates that the
Coulomb interaction plays a key role in the description of magnetism of the
multi-orbital system CuFeAs, especially at the intermediate As height .

Moreover, the antiferromagnetic state with nonvanishing total magnetic moments
is found to be energetically stablized when Cu vacancies are present, which
explains the spontaneous magnetization in nonstoichiometric CuFeAs Thakur  _et
al._ (2014); Qian  _et al._ (2015). Considering the weak ferrimagnetism,
CuFeAs should be susceptible to external magnetic field. It is the reason why
the magnetization exhibits the AFM-like behavior with the decrease of
temperature at low magnetic fields, but shows the FM-like behavior at magnetic
field of $>$ 500 Oe due to the field-induced ferromagnetic component Li  _et
al._ (2018); Thakur  _et al._ (2014); Qian  _et al._ (2015).

Finally, we found that CuFeAs is an unique compound in the family of iron
pnictides which may be used to unveil the origin of weak magnetism present
commonly in the pnictides. In contrast to the CAFM state at lower As height
and the FM state at higher As height, both of which can be understood from the
itinerant electron picture Mazin  _et al._ (2008); Dong  _et al._ (2008);
Mazin (2010) and the localized spin scenario Yildirim (2008); Si and Abrahams
(2008), the BAFM phase which can account for various experimental observations
in CuFeAs Thakur  _et al._ (2014); Qian  _et al._ (2015); Zou  _et al._ (2017)
can only be explained by breakdown of the Hund’s rule. The antiparallel
arrangement of magnetic moments in different orbitals on each iron atom has
been proposed as the possible origin for weak magnetism in LaFeAsO Cricchio
_et al._ (2010); Bascones  _et al._ (2010), which unfortunately was not
prevalent since it was just an alternative theory to the widely accepted ones
based on itinerant or localized scenario Mazin  _et al._ (2008); Dong  _et
al._ (2008); Mazin (2010); Yildirim (2008); Si and Abrahams (2008). However,
CuFeAs may be the first counterexample which casts doubts on the
applicabilities of the well accepted theories and supports breakdown of the
Hund’s rule as a unified picture for the weak magnetism observed
experimentally in iron pnictides. It should be noted that violation of the
first Hund’s rule appears in the presence of two partly filled shells like
cerium Morgan and Kutzelnigg (1993), while there is only one partly filled
shell in CuFeAs. Moreover, Hund’s coupling is always believed to dominate the
correlated metallic behavior in iron pnictides Haule and Kotliar (2009); Yin
_et al._ (2011); Lanatà  _et al._ (2013); Georges  _et al._ (2013) as was
frequently pointed out by dynamical mean field theory (DMFT) Georges  _et al._
(1996) or LDA+DMFT Kotliar  _et al._ (2006) studies where nonlocal
correlations are totally ignored. If breakdown of the Hund’s rule were
dominant, the intersite interorbital hybridizations and multipole interactions
become important, which brings the concept of Hund’s metal into question and
requires further investigations beyond local approximation and Hubbard
interactions. And breakdown of the Hund’s rule may provide a new route to form
singlet cooper pairs locally Hoshino and Werner (2017).

Note that the CuFeAs is unique among iron pnictides due to the fact that it
possesses the highest arsenic height in comparison to other iron arsenic
compounds, for example; $h_{\mathrm{As}}\sim 1.51$ Å for LiFeAs Pitcher  _et
al._ (2008), $1.31$ Å for LaFeAsO de la Cruz  _et al._ (2008), $1.35$ Å for
BaFe2As2 Huang  _et al._ (2008). And it is even higher than that of Fe1.01Se
where $h_{\mathrm{Se}}\sim 1.47$ Å McQueen  _et al._ (2009). The height is
comparable to $h_{\mathrm{Te}}\sim 1.75$ Å of Fe1.068Te which also shows
bicollinear antiferromagnetic order Li  _et al._ (2009b). However, it is
smaller than $h_{\mathrm{Sb}}\sim 1.84$ Å of CuFeSb where ferromagnetism is
observed Qian  _et al._ (2012); Sirohi  _et al._ (2016).

##  V conclusion

In conclusion, we have investigated the magnetism of CuFeAs by applying DFT
calculations. It is found that breakdown of Hund’s rule occurs and is
responsible for the exotic BAFM state in CuFeAs at the height of As atom of
1.612 Å $<h_{\text{As}}<$ 1.84 Å. The novel phase intersects between a CAFM
state at $h<$ 1.612 Å and an FM state at $h>$ 1.84 Å. The presence of Cu
vacancy favors the BAFM state and induce weak ferrimagnetism due to the
symmetry breaking between magnetic sublattices. The interaction is
indispensible to correctly capture the ground state of CuFeAs. Our results can
be applied to fully understand experimental observations and have important
implication that breakdown of Hund’s rule may be a unified theory for weak
magnetism in iron pnictides.

###### Acknowledgements.

This work is financially supported by the National Natural Science Foundation
of China (Grant Nos. 11774258, 12004283) and Postgraduate Education Reform
Project of Tongji University (Grant No. GH1905). Z. Y. Song acknowledges the
financial support by China Postdoctoral Science Foundation (Grant No.
2019M651563).

##  Appendix A  Crystal structure of CuFeAs and antiferromagnetic structures
studied in the work

Figure 5: (color online) (a) Crystal structure of CuFeAs used in this work,
where the blue, brown and green atoms denote Cu, Fe and As, respectively. The
topview cartoons for Néel antiferromagnetic state (NAFM) (b), collinear
antiferromagnetic state (CAFM) (c), and bicollinear antiferromagnetic state
(BAFM) (d), where the red and black arrows represent two different spin
species on Fe atoms which form a square lattice.

## References

  * Kamihara  _et al._ (2008) Yoichi Kamihara, Takumi Watanabe, Masahiro Hirano, and Hideo Hosono, “Iron-Based Layered Superconductor LaO1-xFxFeAs ($x$ = 0.05-0.12) with Tc=26 K,” J. Am. Chem. Soc. 130, 3296 (2008).
  * Rotter  _et al._ (2008) Marianne Rotter, Marcus Tegel, and Dirk Johrendt, “Superconductivity at 38 K in the Iron Arsenide $({\mathrm{Ba}}_{1-x}{\mathrm{K}}_{x}){\mathrm{Fe}}_{2}{\mathrm{As}}_{2}$,” Phys. Rev. Lett. 101, 107006 (2008).
  * Zhi-An  _et al._ (2008) Ren Zhi-An, Lu Wei, Yang Jie, Yi Wei, Shen Xiao-Li, Li Zheng-Cai, Che Guang-Can, Dong Xiao-Li, Sun Li-Ling, Zhou Fang, and Zhao Zhong-Xian, “Superconductivity at 55K in iron-based F-doped layered quaternary compound Sm[O1-xFx]FeAs,” Chin. Phys. Lett. 25, 2215 (2008).
  * Stewart (2011) G. R. Stewart, “Superconductivity in iron compounds,” Rev. Mod. Phys. 83, 1589 (2011).
  * Keimer  _et al._ (2015) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, “From quantum matter to high-temperature superconductivity in copper oxides,” Nature 518, 179 (2015).
  * de la Cruz  _et al._ (2008) Clarina de la Cruz, Q. Huang, J. W. Lynn, Jiying Li, W. Ratcliff, II, J. L. Zarestky, H. A. Mook, G. F. Chen, J. L. Luo, N. L. Wang, and Pengcheng Dai, “Magnetic order close to superconductivity in the iron-based layered LaO1-xFxFeAs systems,” Nature 453, 899 (2008).
  * Basov and Chubukov (2011) D. N. Basov and Andrey V. Chubukov, “Manifesto for a higher T-c,” Nature Phys 7, 272 (2011).
  * Scalapino (2012) D. J. Scalapino, “A common thread: The pairing interaction for unconventional superconductors,” Rev. Mod. Phys. 84, 1383 (2012).
  * de’ Medici  _et al._ (2014) Luca de’ Medici, Gianluca Giovannetti, and Massimo Capone, “Selective Mott Physics as a Key to Iron Superconductors,” Phys. Rev. Lett. 112, 177001 (2014).
  * Boeri  _et al._ (2008) L. Boeri, O. V. Dolgov, and A. A. Golubov, “Is ${\mathrm{LaFeAsO}}_{1-x}{\mathrm{F}}_{x}$ an Electron-Phonon Superconductor?” Phys. Rev. Lett. 101, 026403 (2008).
  * Mazin  _et al._ (2008) I. I. Mazin, D. J. Singh, M. D. Johannes, and M. H. Du, “Unconventional Superconductivity with a Sign Reversal in the Order Parameter of ${\mathrm{LaFeAsO}}_{1-x}{\mathrm{F}}_{x}$,” Phys. Rev. Lett. 101, 057003 (2008).
  * Dong  _et al._ (2008) J. Dong, H. J. Zhang, G. Xu, Z. Li, G. Li, W. Z. Hu, D. Wu, G. F. Chen, X. Dai, J. L. Luo, Z. Fang, and N. L. Wang, “Competing orders and spin-density-wave instability in La(O1-xFx)FeAs,” Europhys. Lett. 83, 27006 (2008).
  * Mazin (2010) Igor I. Mazin, “Superconductivity gets an iron boost,” Nature 464, 183 (2010).
  * Yildirim (2008) T. Yildirim, “Origin of the 150-K Anomaly in LaFeAsO: Competing Antiferromagnetic Interactions, Frustration, and a Structural Phase Transition,” Phys. Rev. Lett. 101, 057010 (2008).
  * Si and Abrahams (2008) Qimiao Si and Elihu Abrahams, “Strong Correlations and Magnetic Frustration in the High ${T}_{c}$ Iron Pnictides,” Phys. Rev. Lett. 101, 076401 (2008).
  * Johannes and Mazin (2009) M. D. Johannes and I. I. Mazin, “Microscopic origin of magnetism and magnetic interactions in ferropnictides,” Phys. Rev. B 79, 220510 (2009).
  * You  _et al._ (2011) Yi-Zhuang You, Fan Yang, Su-Peng Kou, and Zheng-Yu Weng, “Magnetic and superconducting instabilities in a hybrid model of itinerant/localized electrons for iron pnictides,” Phys. Rev. B 84, 054527 (2011).
  * Haule and Kotliar (2009) K. Haule and G. Kotliar, “Coherence-incoherence crossover in the normal state of iron oxypnictides and importance of Hund’s rule coupling,” New J. Phys. 11, 025021 (2009).
  * Yin  _et al._ (2011) Z. P. Yin, K. Haule, and G. Kotliar, “Magnetism and charge dynamics in iron pnictides,” Nature Phys 7, 294 (2011).
  * Lanatà  _et al._ (2013) Nicola Lanatà, Hugo U. R. Strand, Gianluca Giovannetti, Bo Hellsing, Luca de’ Medici, and Massimo Capone, “Orbital selectivity in Hund’s metals: The iron chalcogenides,” Phys. Rev. B 87, 045122 (2013).
  * Georges  _et al._ (2013) Antoine Georges, Luca de’ Medici, and Jernej Mravlje, “Strong Correlations from Hund’s Coupling,” Annu. Rev. Condens. Matter Phys. 4, 137 (2013).
  * Thakur  _et al._ (2014) Gohil S. Thakur, Zeba Haque, L. C. Gupta, and A. K. Ganguli, “CuFeAs: A New Member in the 111-Family of Iron-Pnictides,” J. Phys. Soc. Jpn. 83, 054706 (2014).
  * Li  _et al._ (2018) Kunkun Li, Duanduan Yuan, Jiangang Guo, and Xiaolong Chen, “Observation of direct evolution from antiferromagnetism to superconductivity in $\mathrm{C}{\mathrm{u}}_{1-x}\mathrm{L}{\mathrm{i}}_{x}\mathrm{FeAs}$ ($0\leq x\leq 1.0$),” Phys. Rev. B 97, 134503 (2018).
  * Zou  _et al._ (2017) T. Zou, C. C. Lee, W. Tian, H. B. Cao, M. Zhu, B. Qian, C. R. dela Cruz, W. Ku, Z. Q. Mao, and X. Ke, “$G$-type magnetic order in ferropnictide $\mathrm{C}{\mathrm{u}}_{x}\mathrm{F}{\mathrm{e}}_{1-y}\mathrm{As}$ induced by hole doping on As sites,” Phys. Rev. B 95, 054414 (2017).
  * Kamusella  _et al._ (2017) Sirko Kamusella, Hans-Henning Klauss, Gohil S. Thakur, Zeba Haque, Laxmi C. Gupta, Ashok K. Ganguli, Inga Kraft, Ulrich Burkhardt, Helge Rosner, Hubertus Luetkens, Jeffrey W. Lynn, and Yang Zhao, “Magnetism and site exchange in CuFeAs and CuFeSb: A microscopic and theoretical investigation,” Phys. Rev. B 95, 094415 (2017).
  * Qian  _et al._ (2015) Bin Qian, Jin Hu, Jinyu Liu, Zhida Han, Ping Zhang, Lei Guo, Xuefan Jiang, T. Zou, M. Zhu, C. R. Dela Cruz, X. Ke, and Z. Q. Mao, “Weak ferromagnetism of $\mathrm{Cu}{}_{x}\mathrm{Fe}{}_{1+y}\text{As}$ and its evolution with Co doping,” Phys. Rev. B 91, 014504 (2015).
  * Lv (2009) B. Lv, Ph.D. thesis, University of Houston (2009).
  * Wang  _et al._ (2016) Guangtao Wang, Xianbiao Shi, and Dongyang Wang, “Pnictide-height dependent ferromagnetism in CuFeAs and CuFeSb,”  J. Alloys Compd. 686, 38 (2016).
  * (29) .
  * Tapp  _et al._ (2008) Joshua H. Tapp, Zhongjia Tang, Bing Lv, Kalyan Sasmal, Bernd Lorenz, Paul C. W. Chu, and Arnold M. Guloy, “LiFeAs: An intrinsic FeAs-based superconductor with ${T}_{c}=18\text{ }\text{K}$,” Phys. Rev. B 78, 060505 (2008).
  * Li  _et al._ (2009a) Shiliang Li, Clarina de la Cruz, Q. Huang, G. F. Chen, T.-L. Xia, J. L. Luo, N. L. Wang, and Pengcheng Dai, “Structural and magnetic phase transitions in ${\text{Na}}_{1-\delta}\text{FeAs}$,” Phys. Rev. B 80, 020504 (2009a).
  * Schwarz  _et al._ (2002) K. Schwarz, P. Blaha, and G.K.H. Madsen, “Electronic structure calculations of solids using the WIEN2k package for material sciences,” Comput. Phys. Commun. 147, 71 (2002).
  * Perdew  _et al._ (1996) John P. Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized Gradient Approximation Made Simple,” Phys. Rev. Lett. 77, 3865 (1996).
  * Yin  _et al._ (2010) Wei-Guo Yin, Chi-Cheng Lee, and Wei Ku, “Unified Picture for Magnetic Correlations in Iron-Based Superconductors,” Phys. Rev. Lett. 105, 107004 (2010).
  * (35) .
  * Liechtenstein  _et al._ (1995) A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, “Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators,” Phys. Rev. B 52, R5467 (1995).
  * Czyżyk and Sawatzky (1994) M. T. Czyżyk and G. A. Sawatzky, “Local-density functional and on-site correlations: The electronic structure of ${\mathrm{La}}_{2}$${\mathrm{CuO}}_{4}$ and ${\mathrm{LaCuO}}_{3}$,” Phys. Rev. B 49, 14211 (1994).
  * Miyake  _et al._ (2010) Takashi Miyake, Kazuma Nakamura, Ryotaro Arita, and Masatoshi Imada, “Comparison of Ab initio Low-Energy Models for LaFePO, LaFeAsO, BaFe2As2, LiFeAs, FeSe, and FeTe: Electron Correlation and Covalency,” J. Phys. Soc. Jpn. 79, 044705 (2010).
  * (39) .
  * Yildirim (2009) T. Yildirim, “Strong Coupling of the Fe-Spin State and the As-As Hybridization in Iron-Pnictide Superconductors from First-Principle Calculations,” Phys. Rev. Lett. 102, 037003 (2009).
  * Yin  _et al._ (2008) Z. P. Yin, S. Lebègue, M. J. Han, B. P. Neal, S. Y. Savrasov, and W. E. Pickett, “Electron-Hole Symmetry and Magnetic Coupling in Antiferromagnetic LaFeAsO,” Phys. Rev. Lett. 101, 047001 (2008).
  * Moon and Choi (2010) Chang-Youn Moon and Hyoung Joon Choi, “Chalcogen-Height Dependent Magnetic Interactions and Magnetic Order Switching in ${\mathrm{FeSe}}_{x}{\mathrm{Te}}_{1-x}$,” Phys. Rev. Lett. 104, 057003 (2010).
  * Qian  _et al._ (2012) B. Qian, J. Lee, J. Hu, G. C. Wang, P. Kumar, M. H. Fang, T. J. Liu, D. Fobes, H. Pham, L. Spinu, X. S. Wu, M. Green, S. H. Lee, and Z. Q. Mao, “Ferromagnetism in CuFeSb: Evidence of competing magnetic interactions in iron-based superconductors,” Phys. Rev. B 85, 144427 (2012).
  * Sirohi  _et al._ (2016) Anshu Sirohi, Chandan K. Singh, Gohil S. Thakur, Preetha Saha, Sirshendu Gayen, Abhishek Gaurav, Shubhra Jyotsna, Zeba Haque, L. C. Gupta, Mukul Kabir, Ashok K. Ganguli, and Goutam Sheet, “High spin polarization and the origin of unique ferromagnetic ground state in CuFeSb,” Appl. Phys. Lett. 108, 242411 (2016).
  * Graser  _et al._ (2009) S. Graser, T. A. Maier, P. J. Hirschfeld, and D J Scalapino, “Near-degeneracy of several pairing channels in multiorbital models for the Fe pnictides,” New J. Phys. 11, 025016 (2009).
  * Ding  _et al._ (2013) Ming-Cui Ding, Hai-Qing Lin, and Yu-Zhong Zhang, “Hidden $(\pi,0)$ instability as an itinerant origin of bicollinear antiferromagnetism in Fe1+xTe,” Phys. Rev. B 87, 125129 (2013).
  * Cricchio  _et al._ (2010) Francesco Cricchio, Oscar Grånäs, and Lars Nordström, “Low spin moment due to hidden multipole order from spin-orbital ordering in LaFeAsO,” Phys. Rev. B 81, 140403 (2010).
  * Bascones  _et al._ (2010) E. Bascones, M. J. Calderón, and B. Valenzuela, “Low Magnetization and Anisotropy in the Antiferromagnetic State of Undoped Iron Pnictides,” Phys. Rev. Lett. 104, 227201 (2010).
  * Liu  _et al._ (2012) Chen Liu, Dao-Xin Yao, and Anders W. Sandvik, “Two-orbital quantum spin model of magnetism in the iron pnictides,” Phys. Rev. B 85, 094410 (2012).
  * Mostofi  _et al._ (2008) Arash A. Mostofi, Jonathan R. Yates, Young-Su Lee, Ivo Souza, David Vanderbilt, and Nicola Marzari, “wannier90: A tool for obtaining maximally-localised Wannier functions,”  Comput. Phys. Commun. 178, 685 (2008).
  * Morgan and Kutzelnigg (1993) John D. Morgan and Werner Kutzelnigg, “Hund’s rules, the alternating rule, and symmetry holes,” J. Phys. Chem. 97, 2425 (1993).
  * Georges  _et al._ (1996) Antoine Georges, Gabriel Kotliar, Werner Krauth, and Marcelo J. Rozenberg, “Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions,” Rev. Mod. Phys. 68, 13 (1996).
  * Kotliar  _et al._ (2006) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, “Electronic structure calculations with dynamical mean-field theory,” Rev. Mod. Phys. 78, 865 (2006).
  * Hoshino and Werner (2017) Shintaro Hoshino and Philipp Werner, “Spontaneous Orbital-Selective Mott Transitions and the Jahn-Teller Metal of ${A}_{3}{\mathrm{C}}_{60}$,” Phys. Rev. Lett. 118, 177002 (2017).
  * Pitcher  _et al._ (2008) Michael J. Pitcher, Dinah R. Parker, Paul Adamson, Sebastian J. C. Herkelrath, Andrew T. Boothroyd, Richard M. Ibberson, Michela Brunelli, and Simon J. Clarke, “Structure and superconductivity of LiFeAs,” Chem. Commun. 45, 5918 (2008).
  * Huang  _et al._ (2008) Q. Huang, Y. Qiu, Wei Bao, M. A. Green, J. W. Lynn, Y. C. Gasparovic, T. Wu, G. Wu, and X. H. Chen, “Neutron-Diffraction Measurements of Magnetic Order and a Structural Transition in the Parent ${\mathrm{BaFe}}_{2}{\mathrm{As}}_{2}$ Compound of FeAs-Based High-Temperature Superconductors,” Phys. Rev. Lett. 101, 257003 (2008).
  * McQueen  _et al._ (2009) T. M. McQueen, Q. Huang, V. Ksenofontov, C. Felser, Q. Xu, H. Zandbergen, Y. S. Hor, J. Allred, A. J. Williams, D. Qu, J. Checkelsky, N. P. Ong, and R. J. Cava, “Extreme sensitivity of superconductivity to stoichiometry in ${\text{Fe}}_{1+\delta}\text{Se}$,” Phys. Rev. B 79, 014522 (2009).
  * Li  _et al._ (2009b) Shiliang Li, Clarina de la Cruz, Q. Huang, Y. Chen, J. W. Lynn, Jiangping Hu, Yi-Lin Huang, Fong-Chi Hsu, Kuo-Wei Yeh, Maw-Kuen Wu, and Pengcheng Dai, “First-order magnetic and structural phase transitions in ${\text{Fe}}_{1+y}{\text{Se}}_{x}{\text{Te}}_{1-x}$,” Phys. Rev. B 79, 054503 (2009b).


# Transformation-Dependent Performance-Enhancement of Digital Annealer for
3-SAT

Christian Münch  Fujitsu Services GmbH, Mies-van-der-Rohe-Straße 8, 80807
Munich, Germany  Fritz Schinkel  Fujitsu Services GmbH, Mies-van-der-Rohe-
Straße 8, 80807 Munich, Germany  Sebastian Zielinski  Institute for
Informatics, LMU Munich, Oettingenstraße 67, 80538 Munich, Germany  Stefan
Walter  Fujitsu Services GmbH, Mies-van-der-Rohe-Straße 8, 80807 Munich,
Germany

###### Abstract

Quadratic Unconstrained Binary Optimization (QUBO) problems are NP-hard
problems and many real-world problems can be formulated as QUBO. Currently
there are no algorithms known that can solve arbitrary instances of NP-hard
problems efficiently. Therefore special-purpose hardware such as Digital
Annealer, other Ising machines, as well as quantum annealers might lead to
benefits in solving such problems. We study a particularly hard class of
problems which can be formulated as QUBOs, namely Boolean satisfiability (SAT)
problems, and specifically 3-SAT. One intriguing aspect about 3-SAT problems
is that there are different transformations from 3-SAT to QUBO. We study the
transformations’ influence on the problem solution, using Digital Annealer as
a special-purpose solver. Besides well-known transformations we investigate a
novel in this context not yet discussed transformation, using less auxiliary
variables and leading to very good performance. Using exact diagonalization,
we explain the differences in performance originating from the different
transformations. We envision that this knowledge allows for specifically
engineering transformations that improve a solvers capacity to find high
quality solutions. Furthermore, we show that the Digital Annealer outperforms
a quantum annealer in solving hard 3-SAT instances.

QUBO, Ising, 3-SAT, Annealing, Quantum-Inspired, Optimization

##  I Introduction

Combinatorial optimization deals with discrete decisions and finding the best
combination in a typically huge search space [1, 2]. Such problems are
ubiquitous. Academic problems such as MaxCut, set partitioning, or graph
coloring, as well as many practically relevant problems like task scheduling,
vehicle routing, or knapsack are prime examples of combinatorial optimization
problems. Combinatorial optimization problems have the unpleasant
characteristic in common that the search space grows rapidly (i.e.
exponentially or factorially) with the number of decision variables (i.e.,
possible yes-no decisions). This is the reason why an exhaustive search
approach is very inefficient with increasing problem size. Therefore, many
heuristics, such as simulated annealing [3], tabu search [4], and others are
employed to solve such problems.

A promising framework to solve combinatorial optimization problems is the
Quadratic Unconstrained Binary Optimization (QUBO) approach [5, 6, 7]. QUBO
problems deal with minimizing a polynomial of degree two with 0 and 1 valued
variables. Since there is a correspondence of the QUBO formulation to the
well-known and well-studied Ising spin glass problem [8, 9, 10], the QUBO
formulation gained popularity with the availability of Quantum Annealing
hardware [11, 12, 13, 14, 15, 16] which was quickly followed by other types of
Ising machines, such as Fujitsu’s Digital Annealer [17, 18, 19], Toshiba’s
Simulated Bifurcation machine [20], and others [21].

Besides the aforementioned combinatorial optimization problems, Boolean
satisfiability (SAT) problems (especially 3-SAT) problems are an important
problem class, e.g., to study computational complexity and for instance
showing NP-completeness [22, 23]. SAT problems being constraint satisfaction
problems, are of great practical relevance, with applications mainly in
software verification [26, 27] and hardware verification [28]. There are
various algorithmic approaches on classical CPU and GPU hardware to solve
large and complex SAT problems [29]. However, up to now there is no algorithm
known which is able to solve arbitrary SAT/3-SAT instances in (worst-case)
polynomial time.

3-SAT problems, being combinatorial problems, can also be formulated as QUBO
and novel compute techniques which could provide an advantage over traditional
technologies can be used. In fact, there are studies on (Max-)3-SAT [31, 32,
34, 33], 3XOR-SAT [24], and other variants of satisfiability problems using
quantum annealers, Digital Annealers, and other promising novel compute
architectures which investigate the solvers’ abilities.

In this paper, we use Fujitsu’s Digital Annealer to solve 3-SAT problems and
investigate the performance of different 3-SAT to QUBO transformations. We
explain the differences using small problem instances which allow for exact
diagonalization, and show that extrapolating the findings to larger instances
can be done. We also show that Digital Annealer outperforms a quantum annealer
on solving hard 3-SAT instances.

The paper is structured as follows. In Sec. II we formally introduce the 3-SAT
problem, review existing 3-SAT to QUBO transformations, introduce one novel,
in this context not yet discussed, transformation, and briefly present the
main features of Fujitsu’s Digital Annealer. We continue in Sec. III with our
results obtained using the Digital Annealer. In Sec. IV we present the
analysis of the results. Finally, we conclude and give an outlook in Sec. V.

##  II Background

###  II.1 The 3-SAT Problem and its Relation to QUBO

In general, Boolean satisfiability problems deal with the assignment of
Boolean variables $x_{0},x_{1},...,x_{n-1}$ with $x_{i}\in\\{0,1\\}$ for
$n\in\mathbb{N}$ and $i=0,\ldots,n-1$ such that a Boolean formula $\Phi$ is
satisfied, i.e., $\Phi(x_{0},x_{1},...,x_{n-1})=1$. A Boolean formula $\Phi$
consists of literals (i.e., Boolean variables) and clauses. Literals can be
positive (i.e., they are represented by the variable itself) or negative
(i.e., they are represented by the negated ${\rm{NOT}}$ ($\neg$)) variable.
Literals are connected by the logical ${\rm{OR}}$ ($\vee$) operator and
clauses are joined by the logical ${\rm{AND}}$ ($\wedge$) operator. By
definition a Boolean formula is in conjunctive normal form (CNF) if it is a
conjunction of clauses. Here, we focus on 3-SAT problems in CNF, i.e., clauses
contain exactly three variables. An example with $m=3$ clauses and $n=3$
variables is for instance

$\displaystyle\Phi_{0}=$ $\displaystyle(x_{1}\vee x_{2}\vee\neg x_{3})$
$\displaystyle\wedge$ $\displaystyle(\neg x_{1}\vee x_{2}\vee x_{3})$
$\displaystyle\wedge$ $\displaystyle(\neg x_{1}\vee x_{2}\vee x_{3})$
$\displaystyle\wedge$ $\displaystyle(x_{1}\vee\neg x_{4}\vee x_{5})\,,$ (1)

with the assignment $x_{1}=0,x_{2}=1,x_{3}=1,x_{4}=1,x_{5}=1$ satisfying
$\Phi_{0}$. 3-SAT problems have been proven to be NP-complete [22, 23] which
makes certain instances extremely hard to solve with increasing problem size.
For 3-SAT, in Ref. [25] it was shown that the ratio of clauses to literals
$m/n$ is a good indicator for the hardness of the particular instance. In
particular, 3-SAT instances generated by randomly drawing literals from a
uniform distribution are potentially hard for $m/n\approx 4.24$ [25]. For a
thorough introduction to satisfiability problems we refer to Ref. [30].

As mentioned in Sec. I, combinatorial optimization problems in the form of
QUBOs could potentially be solved efficiently on quantum annealers, gate-based
quantum computers (via the variational Quantum Approximate Optimization
Algorithm), and on other types of Ising machines such as Fujitsu’s Digital
Annealer. The general form of a QUBO is

$\displaystyle
H(\vec{x})=\sum_{i=0}^{n-1}Q_{ii}x_{i}+\sum_{i=0}^{n-2}\sum_{j=i+1}^{n-1}Q_{ij}x_{i}x_{j}\,,$
(2)

where $n\in\mathbb{N},n\geq 2$ defines the problem size and
$Q_{ij}\in\mathbb{R}$ for $i,j\in\\{0,\ldots,n-1\\}$ are the entries of an
$n\times n$ upper triangular matrix ${\bf{Q}}$. Moreover, $\vec{x}$ is a
vector of binary variables $x_{i}\in\\{0,1\\}$ for $i=0,\ldots,n-1$. The QUBO
optimization problem is, to find a vector of binary values for the decision
variables $x_{i}$ such that $\vec{x}$ minimizes $H(\vec{x})$. A mapping from
binary variables $x_{i}$ to spin variables $s_{i}\in\\{-1,1\\}$ via
$s_{i}\mapsto x_{i}=(s_{i}+1)/2$ for $i=0,\ldots,n-1$, transforms the QUBO to
the well-known Ising spin glass problem. Therefore, finding the ground state
of the Ising Hamiltonian
$\tilde{H}(\vec{s})=\sum_{i=0}^{n-1}h_{i}s_{i}+\sum_{i=0}^{n-2}\sum_{j=i+1}^{n-1}J_{ij}s_{i}s_{j}$
(the local fields $h_{i}$ and the couplings $J_{ij}$ define the Ising problem)
is equivalent of solving the QUBO problem [9, 10].

The transformation from 3-SAT to QUBO form has already been studied in
literature. There are various different transformations available: Choi has
introduced a reformulation of 3-SAT instances to instances of a Maximum Weight
Independent Set problem which can then easily be expressed in QUBO form [31].
Chancellor et al. have shown how to express 3-SAT instances as Ising spin
glass models [32]. In Ref. [33], Nüßlein et al. have proposed a method making
use of patterns of the 3-SAT problem to algorithmically construct QUBOs from
3-SAT instances. The approach in Ref. [33] has been further extended to create
thousands of new 3SAT-to-QUBO transformations algorithmically (i.e. by an
automatic procedure, not manually) by Zielinski et al. in Ref. [35]. In this
work, we introduce another, previously not discussed transformation which we
call CountTrue.

For the remainder of the work, we focus on Chancellor, pattern QUBO
transformations, and CountTrue transformations because they show the best
results in terms of solution quality. In the next section we briefly review
the Chancellor, the algorithmic method, to produce thousands of 3SAT-to-QUBO
transformations automatically, and introduce the CountTrue transformation.

###  II.2 Chancellor Transformation

Chancellor’s transformation [32] starts with the binary representation of a
3-SAT formulas’ clause and constructs a sum of parameter dependent QUBOs
representing the satisfiability problem for the formula as a QUBO minimization
problem. In particular, given a binary representation of a clause
$I(\vec{l},\vec{x})=-(l_{1}x_{1}\vee l_{2}x_{2}\vee l_{3}x_{3})$ with
$x_{1},x_{2},x_{3}\in\\{0,1\\}$ and $l_{1},l_{2},l_{3}\in\\{\neg\neg\
,\neg\\}$ expressing the not-negation or negation of the literals, the clause
is then transferred to Ising form. Specifically, for $i=1,2,3$ one transforms
$x_{i}\mapsto s_{i}=\left(2x_{i}-1\right)\in\\{-1,1\\}$ and

$\displaystyle l_{i}\mapsto c_{i}=\begin{cases}1,&\text{if }l_{i}=\neg\neg\\\
-1,&\text{if }l_{i}=\neg\end{cases}\,.$ (3)

This allows to formulate the clause $I(\vec{l},\vec{x})$ in cubic form

$\displaystyle I(\vec{l},\vec{x})\mapsto$
$\displaystyle\frac{1}{8}\tilde{I}_{\rm{clause}}(\vec{c},\vec{s})$
$\displaystyle=$
$\displaystyle-\frac{7}{8}-\frac{1}{8}\left(c_{1}s_{1}+c_{2}s_{2}+c_{3}s_{3}\right)$
$\displaystyle+\frac{1}{8}\left(c_{1}c_{2}s_{1}s_{2}+c_{1}c_{3}s_{1}s_{3}+c_{2}c_{3}s_{2}s_{3}\right)$
$\displaystyle-\frac{1}{8}c_{1}c_{2}c_{3}s_{1}s_{2}s_{3}\,.$ (4)

Equation (II.2) has minimal energy for each satisfying assignment and higher
energy for the only non-satisfying assignment. An auxiliary variable
$s_{a}\in\\{-1,1\\}$ is introduced to reduce the cubic term to quadratic form
and obtain an Ising representation of the clause. In contrast to the proposed
replacement of the cubic term $-c_{1}c_{2}c_{3}s_{1}s_{2}s_{3}$ in Ref. [32]
we choose the replacement

$\displaystyle-c_{1}c_{2}c_{3}s_{1}s_{2}s_{3}\mapsto$ $\displaystyle
I_{\rm{cubic}}(\vec{s},s_{a})$ $\displaystyle\coloneqq$ $\displaystyle
J\sum_{i=1}^{3}\sum_{j=1}^{i-1}s_{i}s_{j}-c_{1}c_{2}c_{3}\sum_{i=1}^{3}s_{i}$
$\displaystyle+2J\sum_{i=1}^{3}s_{i}s_{a}-2c_{1}c_{2}c_{3}s_{a}\,,$ (5)

where $J\geq 1$ is a free parameter. Although not explicitly mentioned in Ref.
[35], Eq. (II.2) is the replacement used in their experiments. Note that with
given fixed values for $s_{1}^{*},s_{2}^{*},s_{3}^{*}$ and free variable
$s_{a}$, minimizing $I_{\rm{cubic}}(s_{1}^{*},s_{2}^{*},s_{3}^{*},s_{a})$ over
$s_{a}$ is a parity check for the variables
$c_{1}s_{1}^{*},c_{2}s_{2}^{*},c_{3}s_{3}^{*}$ (the $c_{i}$ are fixed by the
actual clause). As a result of the minimization over $s_{a}$, the minimal
energy is $E_{\rm{min}}<0$ for uneven numbers of ones $(+1)$, i.e.,
$\\{c_{1}s_{1}^{*},c_{2}s_{2}^{*},c_{3}s_{3}^{*}\\}\in\\{\\{-1,-1,1\\},\\{1,1,1\\}\\}$,
and $E_{\rm{min}}+2$ for even numbers of ones, i.e.,
$\\{c_{1}s_{1}^{*},c_{2}s_{2}^{*},c_{3}s_{3}^{*}\\}\in\\{\\{-1,-1,-1\\},\\{-1,1,1\\}\\}$.
The remaining terms in Eq. (II.2) imply a minimal energy
$-7+E_{\rm{min}}+2+6=E_{\rm{min}}+1$ for the non-satisfying choice
$(c_{1}s_{1}^{*},c_{2}s_{2}^{*},c_{3}s_{3}^{*})=(-1,-1,-1)$ and a minimal
energy $-7+E_{\rm{min}}$ for all satisfying choices. This results in the
representation of the clause satisfiability problem as minimization problem of

$\displaystyle I_{\rm{clause}}(\vec{s},s_{a})=$
$\displaystyle-7-\left(c_{1}s_{1}+c_{2}s_{2}+c_{3}s_{3}\right)$
$\displaystyle+\left(c_{1}c_{2}s_{1}s_{2}+c_{1}c_{3}s_{1}s_{3}+c_{2}c_{3}s_{2}s_{3}\right)$
$\displaystyle+I_{\rm{cubic}}(\vec{s},s_{a})\,.$ (6)

Transforming to binary variables $s_{i}\mapsto x_{i}=(s_{i}+1)/2$ for
$i\in\\{1,2,3,a\\}$ we obtain a QUBO representing the clause satisfiability
problem. Adding up the corresponding QUBOs for each clause (with separate
auxiliary variable for each clause) we obtain a $J$-dependent QUBO
representing the satisfiability problem for the formula.

In the experiments of this paper we choose the free parameter as $J=1$ and
$J=5$ and call the resulting transformations ChancellorJ1 and ChancellorJ5,
respectively.

###  II.3 Automated Pattern QUBOs

In Ref. [35] a method for automatically creating new QUBOs, exploiting
properties of Pattern QUBOs is presented. To understand this approach we first
need to observe, that there are only 4 types of clauses that can occur in a
3SAT problem:

  * •

Type 0 := $(a\vee b\vee c)$

  * •

Type 1 := $(a\vee b\vee\neg c)$

  * •

Type 2 := $(a\vee\neg b\vee\neg c)$

  * •

Type 3 := $(\neg a\vee\neg b\vee\neg c).$

These types of clauses are created by reordering a given 3-SAT clause, such
that all negated literals are at the end of the clause. Now, all that is left
to create a (new) 3-SAT to QUBO transformation is to create a QUBO for each of
the 4 types of clauses. That is, a 3-SAT to QUBO transformation in this
approach consists of exactly 4 reusable QUBOs, which are called Pattern QUBOs
- one for each type of clause. The ground states of Pattern QUBOs should
correspond to correct solutions of a clause. As there are 3 bits in a clause,
there are $2^{3}=8$ possible bit strings per clause, of which 7 are satisfying
and only 1 is not satisfying the clause. Therefore, Pattern QUBO that
corresponds to a given clause type should have at least one ground state for
each of the satisfying assignments and the not satisfying assignment should
correspond to an excited state. After finding a Pattern QUBO for each of the
clause types one proceeds to create a QUBO for the whole 3-SAT problem. This
is done, by iterating over all the clauses from the 3-SAT problem and
transforming all the clauses of the problem into QUBOs (using one of the 4
previously created Pattern QUBOs) and then combining all these clause QUBOs
into a single large QUBO by adding them up using the technique of
superimposing. This way, a ground state in the large QUBO directly corresponds
to a solution of the given 3-SAT problem.

This approach begs an immediate question: As each 3-SAT to QUBO transformation
consists of 4 Pattern QUBOs, how can we find these Pattern QUBOs? In Ref. [35]
the authors present a simple algorithmic method that solves this problem. To
be able to create a Pattern QUBO for a given 3-SAT clause that satisfies the
constraint that all satisfying assignments of a 3-SAT clause need to
correspond to at least one ground state in the Pattern QUBO, the Pattern QUBO
needs to be of dimension 4 (i.e., a $4\times 4$ upper triangular matrix. See
Ref. [35] for more details). Since each Pattern QUBO is a $4\times 4$ upper
triangular QUBO matrix, there are only 10 fields of values to be filled in.
This is small enough for an exhaustive search algorithm to find significant
amounts of solutions. Thus, given a set of values to choose from, the
algorithm inserts all possible combinations of the given set of values into
the 10 fields of the $4\times 4$ Pattern QUBO matrix, to find Pattern QUBOs
that correspond to one of the 4 clause types (i.e. that has at least one
ground state for each satisfying assignment in the clause). Note that the
Chancellor transformation can also be created automatically using this
approach. Thus, because ChancellorJ1 and ChancellorJ5 behave very differently
(see Ref. [34]), the performance of QUBO transformations created by this
approach heavily depends on the choice of Pattern QUBOs used for the different
types of clauses.

As QUBOs generated by this approach can lead to better results, when solving
3-SAT problems on quantum annealers (as shown in Ref. [35]), we also used the
best performing QUBO of Ref. [35]. We will call this QUBO the _AlgorihmQUBO_
for the remainder of this paper. The AlgorithmQUBO consists of four Pattern
QUBOs shown in Tab. 1

|  a |  b |  c |  K
---|---|---|---|---
a |  |  1 |  |  -1
b |  |  |  |  -1
c |  |  |  -1 |  1
K |  |  |  |

Type 0 - Pattern QUBO

|  a |  b |  c |  K
---|---|---|---|---
a |  |  |  |  -1
b |  |  |  -1 |  1
c |  |  |  1 |  -1
K |  |  |  |  1

Type 1 - Pattern QUBO

|  a |  b |  c |  K
---|---|---|---|---
a |  |  -1 |  |  1
b |  |  1 |  |  -1
c |  |  |  |  1
K |  |  |  |

Type 2 - Pattern QUBO

|  a |  b |  c |  K
---|---|---|---|---
a |  |  |  |  1
b |  |  |  1 |  -1
c |  |  |  |  -1
K |  |  |  |  1

Type 3 - Pattern QUBO

Table 1:  Pattern QUBOs for the four different types of clauses of the
AlgorithmQUBO [35].

The variables $K$, in all of the Pattern QUBOs in Tab. 1 represent additional,
so-called _ancilla variables_. They have no representation in the
corresponding 3-SAT problem but need to be added to correctly represent a
3-SAT clause as a QUBO optimization problem (see Ref. [35] for further
details). Note, that for each clause a new ancilla variable needs to be added.
Thus if there are $m$ clauses in the 3-SAT problem, this transformation needs
$m$ additional variables.

###  II.4 CountTrue: A 3-SAT Transformation with HOBO to QUBO Reduction

We now introduce a straightforward approach to cast 3-SAT problems into QUBO
form which is based on quadratization [36, 37, 38, 39, 40, 41]. The approach
we follow is to penalize (per clause) variable configurations which do not
satisfy a clause. It can also be viewed as counting the ${\rm{TRUE}}$ values
per clause (since clauses with one, two, or three ${\rm{TRUEs}}$ are
satisfied) and assign energy zero if the number of ${\rm{TRUEs}}$ is larger
than zero and a positive energy value if there is no ${\rm{TRUE}}$ in the
clause. This way, a solution with zero energy satisfies the 3-SAT problem,
whereas a solution with energy larger than zero does not satisfy the problem.
If we take a clause $(l_{1}x_{1}\vee l_{2}x_{2}\vee l_{3}x_{3})$ with
$x_{1},x_{2},x_{3}\in\\{0,1\\}$ and $l_{1},l_{2},l_{3}\in\\{\neg\neg\
,\neg\\}$, we would need to construct a binary polynomial which penalizes
configurations where $l_{1}x_{1}=l_{2}x_{2}=l_{3}x_{3}=0$. In general a
polynomial with the aforementioned properties is the following

$\displaystyle H_{\rm{CT}}(\vec{x})=-\sum_{c=0}^{m-1}$
$\displaystyle\left(\sum_{(x,\chi_{x})\in\mathcal{I}_{c}}\left(\frac{1-\chi_{x}}{2}+\chi_{x}x\right)-1\right)$
(7)
$\displaystyle\cdot\left(\sum_{(x,\chi_{x})\in\mathcal{I}_{c}}\left(\frac{1-\chi_{x}}{2}+\chi_{x}x\right)-2\right)$
$\displaystyle\cdot\left(\sum_{(x,\chi_{x})\in\mathcal{I}_{c}}\left(\frac{1-\chi_{x}}{2}+\chi_{x}x\right)-3\right)\,,$

where $\mathcal{I}_{c}$ is the set of bits $x$ and corresponding signs
$\chi_{x}$ of literals in clause $c$. In particular, $\chi_{x}$ is the sign of
the literal in the clause, i.e., $\chi_{x}=1$ if $x$ is in the clause and
$\chi_{x}=-1$ if $\neg x$ is in the clause. Unfortunately, $H_{\rm{CT}}$ is a
higher order binary polynomial (HOBO) because it contains cubic terms. A
reduction to QUBO is straight forward by replacing a product of two bits by a
new variable (e.g., $y\in\\{0,1\\}$) and coupling the value of the new
variable to the value of the product by a penalty polynomial. The penalty
polynomial then needs to attain the value zero iff the new variable has the
same value as the product of the two bits. For example, in the product of
three variables $x_{1}x_{2}x_{3}$ we replace $x_{2}x_{3}$ by $y$ where $y$ is
a new bit variable. To enforce that $y=1$ iff $x_{1}x_{2}=1$ we need to add an
additional penalty term of the form $x_{1}x_{2}-2x_{1}y-2x_{2}y+3y$. It is
worth pointing out that compared to the other transformations, CountTrue uses
the minimal amount of auxiliary bits since they can often be reused in several
clauses. See also Fig. 1 for an overview of the bits needed for each of the
transformations.

###  II.5 Digital Annealer

Besides Quantum Annealing and variational algorithms on gate-based quantum
computers, Fujitsu’s Digital Annealer Unit (DAU) is a very promising and
powerful technology to solve large and complex combinatorial optimization
problems in form of QUBOs. The DAU is a custom application-specific integrated
circuit (ASIC) hardware architecture realized using conventional CMOS
technology. A Markov Chain Monte Carlo (MCMC) algorithm forms the basis of the
search algorithm. By implementing the MCMC algorithm on hardware, the search
for low energy solutions to an optimization problem can drastically be
enhanced. Specifically, using techniques such as parallel search, a dynamic
offset energy to escape from local minima, and parallel tempering speeds up
the search. With version 2 of the DAU (DA V2) problem instances with up to
8192 fully connected variables and an integer precision of 16 $\rm{bit}$ for
the QUBO values $Q_{ij}$ can be treated. If the integer precision needs to be
increased to 64 $\rm{bit}$, the maximal problem size is reduced to 4096. In
comparison to quantum annealers, no special environment such as cryogenics,
vacuum, etc. is needed. Another limiting factor present in current quantum
annealers is the so-called minor embedding which is absent in Digital Annealer
due to the all-to-all connectivity of decision variables. In addition, problem
sizes which are of practical relevance and go beyond many toy model
considerations when using quantum computing can be solved with Digital
Annealer.

QUBO implementation and addressing the solver are done using Futjitsu’s
proprietary software development kit (SDK) called `dadk` [42]. Further details
on the Digital Annealer, the SDK, solver parameter, and use cases can also be
found in Ref. [42].

##  III Results

In the following, we introduce and characterize the 3-SAT instances and
furthermore present results in terms of percentage of solved 3-SAT instances
and percentage of correct solutions using DA V2.

###  III.1 3-SAT Instances

We use the same 3-SAT instances with $m=46$ clauses and $n=11$ variables as in
Ref. [35] to evaluate the performance of the different transformations. The
instances were randomly generated drawing from a uniform distribution. With
$m/n=46/11\approx 4.18$ the ratio is chosen close the 3-SAT phase transition
which makes them potentially hard to solve [25]. In total we investigate 1000
different 3-SAT instances. As also pointed out in Ref. [35], it proofs
insightful to inspect properties of the QUBO matrix for each of the
transformations, namely the number of bits needed, the number of different
coefficient values of quadratic terms, and the distribution of the values
range of the coefficient values of quadratic terms. We show these properties
in Fig. 1. We see that CountTrue needs the lowest amount of bits which can be
explained by the fact that we reuse auxiliary variables (for some instances
auxiliary variables can be reused more often for other not). From Fig. 1 we
also see that ChancellorJ5 and CountTrue, as well as AlgorithmQUBO and
ChancellorJ1 show similar structure in their QUBO matrices.

Figure 1:  QUBO matrix properties for each of the transformations. The plots
show a boxplot characterizing all 1000 different 3-SAT instances. a) Number of
bits needed to represent the 3-SAT instance as a QUBO. Here, CountTrue varies
in the number of bits which is due to the fact that we reuse auxiliary
variables. b) Number of different quadratic values of the QUBO matrix. c) Size
of the value range of quadratic values, i.e., the absolute value of difference
largest quadratic value minus smallest quadratic value of the QUBO matrix.

Furthermore, using the same instances as in Ref. [35] allows us to compare the
results obtained in Ref. [35] with a D-Wave quantum annealer to the ones we
obtain using Fujitsu’s Digital Annealer Version 2.

###  III.2 Experiment - Settings

We use Digital Annealer V2 to solve each of the 1000 3-SAT instances 100
times. We make use of an automatic scaling which scales the entries of the
QUBO matrix in a way that the largest value of the QUBO matrix is less than
$2^{p-1}$ where $p=64$ is the integer precision of coefficients in the QUBO
matrix. Furthermore, we use an automatic determination of the annealing
temperatures, i.e., the start and end temperature as well as the dynamical
offset energy; details can be found in Ref. [42] and App. A. We use DA V2 in
annealing mode, not in parallel tempering mode. For each transformation we
vary the number of iterations of DA V2 from $10^{3}$ to $10^{8}$. In the case
of the CountTrue transformation, we need to reduce a binary polynomial of
third degree to second degree. This is conveniently done using the
`reduce_higher_degree` method of the `BinPol` class in the `dadk` SDK [42].

###  III.3 Experiment - Results

Figure 2 shows the results of the experiment. In Fig. 2 a) we show the
percentage of solved 3-SAT instances and in Fig. 2 b) we show the percentage
of overall correctly solved 3-SAT instances, i.e., out of $100\times 1000$
experiments for each fixed iteration number and all transformations.

Figure 2:  a) The percentage of solved 3-SAT instances increases with
increasing number of iterations for all transformations. Vertical lines
indicate the minimal number of iterations to reach $100\%$ correct solutions.
b) The percentage of correct solutions also increases with the number of
iterations. Two groups are forming among the transformations, showing similar
behavior, namely, AlgorithmQUBO / ChancellorJ1 and CountTrue / ChancellorJ5.
This grouping is also reflected in the properties of the QUBO matrix, see Fig.
1. CountTrue shows the best performance which we attribute to the smaller
amount of bits needed to represent the 3-SAT instances as QUBO.

We see that with increasing number of iterations, the percentage of solved
3-SAT instances and the percentage of correct solutions also increase.

Furthermore, we see that with all transformations it is possible to solve all
100 different 3-SAT instances [the vertical dashed colored lines in Fig. 2 a)
indicate the minimal iteration number where $100\%$ of the 3-SAT instances are
solved]. However, the probability for finding the correct solution is (in most
cases) still less than $100\%$ at this point, but for some transformations
increases up to $100\%$ at higher iterations, cf. Fig. 2 b).

In addition, we see from Fig. 2 b) that for less iterations AlgorithmQUBO and
ChancellorJ1 lead to a substantial amount of correct solutions compared to
CountTrue and ChancellorJ5. However, there is a crossover: with increasing
number of iterations CountTrue and ChancellorJ5 outperform the other two
transformations [also see inset of Fig. 2 b)]. Actually, CountTrue and
ChancellorJ5 are the only two transformations reaching $100\%$ correct
solutions [indicated by the vertical dashed colored lines in Fig. 2 b)] at
all.

The similar behavior of CountTrue and ChancellorJ5 on the one hand and
AlgorithmQUBO and ChancellorJ1 on the other hand, is also reflected in their
similar properties of the QUBO matrix, see Fig. 1. However, there are
differences between CountTrue and ChancellorJ5: CountTrue shows an earlier
crossover behavior and reaches $100\%$ correct solution before ChancellorJ5.
This makes CountTrue perform better than ChancellorJ5 for large numbers of
iterations. We attribute this fact to the smaller number of bits needed for
CountTrue than for ChancellorJ5.

In Ref. [35] the D-Wave Advantage_system4.1 was used to solve the same problem
instances and the transformations ChancellorJ1 and ChancelorJ5 were also
investigated. Comparing our results which we obtained with DA V2, we see that
already with as few as $10^{4}$ iterations DA V2 achieves better results than
the D-Wave quantum annealer. With increasing number of iterations, DA V2
drastically outperforms the D-Wave quantum annealer.

| Solved 3-SAT instances $[\%]$ | Correct solutions $[\%]$
---|---|---
| D-Wave | DA V2 | D-Wave | DA V2
DA iterations |  | $10^{4}$ | $10^{8}$ |  | $10^{4}$ | $10^{8}$
ChancellorJ1 | 91.0 | 99.3 | 100.0 | 8.71 | 21.27 | 99.7
ChancellorJ5 | 58.0 | 61.1 | 100.0 | 0.48 | 1.16 | 100.0
Table 2:  Comparison between D-Wave quantum annealer and DA V2. For DA V2 we
show the results for $10^{4}$ and $10^{8}$ iterations. Numbers for D-Wave are
taken from Refs. [35] where the results were obtained using default
parameters. Also, from Fig. 2 it is obvious that with increasing iterations,
the solution quality of DA V2 is much better than using a D-Wave quantum
annealer.

The D-Wave results were obtained using the default parameter settings with an
anneal time of $20\,{\rm{\mu s}}$. For DA V2, $10^{4}$ iterations result in
$\approx 16\,{\rm{ms}}$ anneal time on the DA V2 chip. For $10^{6}$ iterations
this time increases to $\approx 142\,{\rm{ms}}$ anneal time and for $10^{8}$
iterations the anneal time on the DA V2 chip is $\approx 12.8\,{\rm{s}}$.
Times are for obtaining all $100$ solutions. Since the anneal time on the
quantum annealer is limited (on the hardware side) by the qubit coherence
time, a comparison on the time to solution level is hard.

##  IV Analysis of Transformations and their Performances

In order to better understand the differences in performance of the
transformations with regards to the solution quality, especially the
percentage of correct solutions, it would be very instructive to know more
about the energy landscape of the QUBO formulation resulting from the
different transformations.

Although the QUBOs are very small ($<100$ bits), exact diagonalization is out
of reach to learn more about eigenvalues and eigenstates. This is due to the
fact that the corresponding Hamiltonian matrix grows exponentially with the
number of bits needed to represent the 3-SAT instance as a QUBO. Therefore, we
study small 3-SAT instances and diagonalize the corresponding Ising
Hamiltonian such that we are able to investigate the eigensystem of the
Hamiltonian matrix. Here, we take ten randomly created 3-SAT instances with
$m=20$ clauses and $n=5$ variables. The properties of the QUBO matrices are
shown in Fig. 3 a)-b). As in the case of the large instances (see Fig. 1),
also for the small instances, we encounter two groups, AlgorithmQUBO /
ChancellorJ1 and CountTrue / ChancellorJ5 which show similar QUBO properties.

In the next sections, we analyze the small instances exactly, explore
correlations to the solution quality, and extrapolate the gained insights to
the larger 3-SAT instances.

###  IV.1 Exact Diagonalization

Diagonalization of the QUBO matrix without constraining to binary variables
would lead to solutions which are not meaningful to the actual problem.
Therefore, we apply the transformation
$F(\vec{x}):=\bigotimes_{i=0}^{n-1}f_{i}(x_{i})$, with

$\displaystyle f_{i}:\\{0,1\\}\rightarrow\mathbb{R}^{2},\quad
u\mapsto\left(\begin{matrix}u\\\ 1-u\end{matrix}\right),$ (8)

to the binary vector $\vec{x}\in\\{0,1\\}^{n}$. We then examine the QUBO
problem in the vector space $\mathbb{R}^{2^{n}}$. In particular, we elevate
each binary variable $x_{i}$ occurring in the QUBO to a matrix by mapping
$x_{i}\mapsto(I_{i}+Z_{i})/2$ where $I$ is the $2\times 2$ identity matrix and
$Z={\rm{diag}}(1,-1)$ is the Pauli-$z$ matrix. Identity matrices have to be
inserted appropriately when transforming. For example, in a QUBO with $n$
binary variables $x_{i}$, $i=0,\ldots,n-1$, linear terms $\sim x_{i}$ are
mapped by

$\displaystyle
x_{i}\mapsto\bigotimes_{0<\alpha<i}I_{\alpha}\otimes\frac{I_{i}+Z_{i}}{2}\otimes\bigotimes_{i<\beta\leq
n-1}I_{\beta}\,,$ (9)

and quadratic terms $\sim x_{i}x_{j}$ are mapped in the following way

$\displaystyle x_{i}x_{j}\mapsto$
$\displaystyle\bigotimes_{0<\alpha<i}I_{\alpha}\otimes\frac{I_{i}+Z_{i}}{2}\otimes$
(10)
$\displaystyle\bigotimes_{i<\beta<j}I_{\beta}\otimes\frac{I_{j}+Z_{j}}{2}\otimes$
$\displaystyle\bigotimes_{j<\gamma\leq n-1}I_{\gamma}\,.$

The total Hamiltonian matrix is then the sum of all transformed QUBO terms and
is by construction diagonal.

For the small instances we construct the QUBO matrix as before, do the mapping
according to Eqs. (9) and (10), and exactly diagonalize the resulting
Hamiltonian matrix. First, we have a closer look at the degeneracy of the
smallest eigenvalue of the spectrum.

Figure 3:  QUBO matrix properties for each of the transformations for small
3-SAT instances $m=20$ and $n=5$ . The plots show a boxplot characterizing all
10 different 3-SAT instances. a) Number of different quadratic values of the
QUBO matrix. b) Size of the value range of quadratic values, i.e., the
absolute value of difference between largest quadratic value and smallest
quadratic value of the QUBO matrix. c) Energy gap between the smallest
eigenvalue and the second smallest eigenvalue. d) Degeneracy of the smallest
eigenvalue (orange) and second smallest eigenvalue (light orange). As in the
case for the other instances (cf. Fig. 1) two groups AlgorithmQUBO /
ChancellorJ1 and CountTrue / ChancellorJ5 can be distinguished.

In Fig. 3 d) we show the degeneracy of the smallest and second smallest
eigenvalue for each transformation. Again, we see two groups, the first group
AlgorithmQUBO / ChancellorJ1 shows a relatively high degeneracy (both for the
smallest and second smallest eigenvalue), in contrast to that the second group
CountTrue / ChancellorJ5 shows a relatively low degeneracy. This grouping also
exists for the energy gap between the smallest and second smallest
eigenvalues, shown in Fig. 3 c). Except for the CountTrue transformation, the
search space is the same.

Since, the degeneracy of the smallest eigenvalue is directly related to the
number of solutions of the problem instance, a high ground state degeneracy
(at constant size of the search space) leads to an increased probability of
finding this minimum energy state. This explains the behavior of the
transformations shown in Fig. 2. At already low number of iterations, the
transformations AlgorithmQUBO / ChancellorJ1 yield quite good solutions, which
is due to the fairly large number of correct solutions (i.e., high degeneracy
of the smallest eigenvalue) with respect to the size of the search space.
Increasing the number of iterations saturates the curves for AlgorithmQUBO and
ChancellorJ1 in Fig. 2 b). For CountTrue / ChancellorJ5, however, it is harder
to find a minimum energy solution with just a few iterations because the
number of solutions is very small compared to the size of the search space.
Increasing the number of iterations increases the probability to find a
minimum energy solutions for CountTrue / ChancellorJ5.

Moreover, the transformations AlgorithmQUBO / ChancellorJ1 also have a high
degeneracy of the second lowest eigenvalue which are of low energy [see the
gap shown in Fig. 3 c)] but which do not satisfy the SAT instances. This fact
might hinder the improvement in correct solutions with increasing iterations,
since many iterations might be spent on flipping back and forth between
configurations corresponding to the second lowest eigenvalue (which have the
same energy value) rather than finding configurations with minimal energy
value (i.e. correct solutions).

In contrast, the CountTrue / ChancellorJ5 transformations show a small number
of correct solutions for small iterations which is due to the low degeneracy
of the lowest eigenvalue. However, they show good convergence towards the
optimal solution with increasing iterations which is due to the low degeneracy
of the second lowest eigenvalue. In particular, for CountTrue / ChancellorJ5,
there are fewer energetically semi-good but incorrect configurations (because
these correspond to the second lowest eigenstate) in which the annealing
process might get stuck.

From Fig. 2 b), we see that CountTrue and ChancellorJ5 do not behave as
similar as AlgorithmQUBO / ChancellorJ1. To be precise, the curve (c.f. Fig. 2
b) black and red lines) for the CountTrue transformation is slightly shifted
to the left in comparison to the curve for the ChancellorJ5 transformation. We
attribute this shift to the smaller search space of the CountTrue
transformation which is due to the smaller number of bits needed.

To summarize, the insights from analyzing small instances exactly can be
extrapolated to the large instances and provide a reasonable ground to explain
differences in the transformations on the solution of 3-SAT instances.

###  IV.2 Hamming Distance

To further strengthen the argument we investigate the difference of the
eigenstate configurations (i.e., solution configurations corresponding to the
eigenstates of minimal eigenvalue) in terms of the Hamming distance between
eigenstate configurations. We show the statistics of the Hamming distance for
the eigenstate configurations corresponding to the smallest eigenvalues, the
eigenstate configurations corresponding to the second smallest eigenvalues,
and the eigenstate configurations corresponding to the smallest and the second
smallest eigenvalues in Fig. 4.

Figure 4:  Average frequency of the Hamming distance among eigenstate
configurations corresponding to the smallest eigenvalues, the eigenstate
configurations corresponding to the second smallest eigenvalues, and the
eigenstate configurations corresponding to the smallest and the second
smallest eigenvalues (left panels to right panels) for each transformation [a)
to d)]. We see a grouping of AlgorithmQUBO / ChancellorJ1 and CountTrue /
ChancellorJ5 explaining their similar influence on the solution. However,
CountTrue and ChancellorJ5 show slightly different statistics which is due to
their different search spaces.

Again, we see a grouping of AlgorithmQUBO / ChancellorJ1 and CountTrue /
ChancellorJ5. This time however with subtle differences between the CountTrue
and ChancellorJ5 transformation.

As a first observation, we see that the average frequency of Hamming distance
between eigenstate configurations of smallest eigenvalues and also of the
second smallest eigenvalues is higher for the ChancellorJ1 transformation
compared to the AlgorithmQUBO transformation. In particular, AlgorithmQUBO and
ChancellorJ1 show a peak at Hamming distance eight between different
configurations corresponding to the second smallest eigenvalue. The peak is
about 1.5 times higher for ChancellorJ1. Extrapolating to the larger 3-SAT
instance and the result presented in Fig. 2, this explains the slightly slower
growth of ChancellorJ1 in Fig. 2 which can then be associated to spending many
iterations on flipping between configurations corresponding to the second
smallest eigenvalue.

Concerning the CountTrue and ChancellorJ5 transformation, we see a similar
degeneracy for smallest and second smallest eigenvalues. However, there are
slight differences in the average Hamming distance frequencies between
configurations corresponding to smallest, second smallest, and smallest to
second smallest eigenvalues. The average frequency for Hamming distance
between configurations corresponding to the smallest eigenvalue is always
lower than one for both transformations. This is due to the fact that for most
formulas there are only a single or very few valid solutions for CountTrue and
ChancellorJ5. This explains the similar behavior of both transformations for
small iteration numbers with a slight advantage for CountTrue probably due to
the smaller search space. Concerning the second smallest eigenvalues and
corresponding configurations, the average frequency for small Hamming
distances is slightly higher for CountTrue than for ChancellorJ5. In
principle, we would therefore expect a slower increase in correct solutions
with respect to iteration numbers for CountTrue. We account the opposite
behavior in Fig. 2 for low iteration numbers to the smaller search space of
CountTrue which seems to compensate the expected slower convergence due to
smaller Hamming distances of configurations corresponding to the second
smallest eigenvalue. At the same time, the degeneracy for the second smallest
eigenvalue is very low both for CountTrue and ChancellorJ5 which also
mitigates the effect of Hamming distances of configurations corresponding to
the second smallest eigenvalue. However, the convergence above $80\%$ solved
instances / correct solutions is indeed slower for CountTrue compared to
ChancellorJ5 as seen in Fig. 2, which indicates that there exist instances for
which the smaller Hamming distance between configurations corresponding to the
second smallest eigenvalue is problematic.

##  V Conclusion

To summarize, we have studied different 3-SAT to QUBO transformations and
their influence on the solution quality when solving 3SAT instances with
Fujitsu’s Digital Annealer. In addition to the previously studied 3-SAT to
QUBO transformations, we have introduced a novel transformation, CountTrue. We
have seen that different transformations have different impact on the solution
quality. We explained the different behavior by means of exact diagonalization
of small 3-SAT instances and extrapolated the insights to larger instances.
This analysis also explains the good performance of our novel CountTrue
transformation. Furthermore, we have compared our results to results obtained
with D-Wave’s quantum annealer Advantage 4.1. This comparison shows a clear
advantage of Digital Annealer over the quantum annealer in terms of solution
quality. It is worth mentioning that the performance ranking of different
transformations is the same on the quantum annealer and on the Digital
Annealer.

With these insights, we plan to investigate tailoring of ground state
degeneracies to enhance the solution quality. In addition, we plan to
investigate larger 3-SAT instances.

Acknowledgements.– We would like to acknowledge stimulating discussions with
Markus Kirsch, and valuable comments by Jan Budich and Andreas Rohnfelder.

## References

  * [1] L. A. Wolsey and G. L. Nemhauser, Integer and Combinatorial Optimization (Wiley, 1999).
  * [2] C. H. Papadimitriou and K. Steiglitz, Combinatorial Optimization: Algorithms and Complexity, Dover Books on Computer Science (Dover Publications, 2013)
  * [3] S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi, Science 220, 671 (1983).
  * [4] F. Glover, Computers and Operations Research 13, 533 (1986).
  * [5] Z. Lü, F. Glover, and J.K. Hao European Journal of Operational Research 207, 1254 (2010).
  * [6] G. Kochenberger, J.K. Hao, F. Glover, M. Lewis, Z. Lu, H. Wang, and Y. Wang, Journal of Combinatorial Optimization 28, 58 (2014).
  * [7] F. Glover and G. Kochenberger, arXiv:1811.11538 (2019).
  * [8] F. Barahona, M. Grötschel, M. Jünger, and GṘeinelt, Operations Research 36, 493 (1988).
  * [9] A. Lucas, Frontiers in Physics 2, 5 (2014).
  * [10] S. Boettcher, Phys. Rev. Research 1, 033142 (2019).
  * [11] M. W. Johnson, P. Bunyk, F. Maibaum, E. Tolkacheva, A. J. Berkley, E. M. Chapple, R. Harris, J. Johansson, T. Lanting, I. Perminov, E. Ladizinsky, T. Oh, and G. Rose, Superconductor Science and Technology 23, 065004 (2010).
  * [12] A. J. Berkley, M. W. Johnson, P. Bunyk, R. Harris, J. Johansson, T. Lanting, E. Ladizinsky, E. Tolkacheva, M. H. S. Amin, and G. Rose, Superconductor Science and Technology 23, 105014 (2010).
  * [13] R Harris, M. W. Johnson, T. Lanting, A. J. Berkley, J. Johansson, P. Bunyk, E. Tolkacheva, E. Ladizinsky, N. Ladizinsky, T. Oh, F. Cioata, I. Perminov, P. Spear, C. Enderud, C. Rich, S. Uchaikin, M. C. Thom, E. M. Chapple, J. Wang, B. Wilson, M. H. S. Amin, N. Dickson, K. Karimi, B. Macready, C. J. S. Truncik, and G. Rose, Phys. Rev. B 82, 024511 (2010).
  * [14] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, Nature 473, 194 (2011).
  * [15] P. I. Bunyk, E. M. Hoskinson, M. W. Johnson, E. Tolkacheva, F. Altomare, A. J. Berkley, R. Harris, J. P. Hilton, T. Lanting, A. J. Przybysz, and J. Whittaker, IEEE Transactions on Applied Superconductivity 24, 1 (2014).
  * [16] K. Boothby, P. Bunyk, J. Raymond, and A. Roy, arXiv:2003.00133 (2020)
  * [17] S. Matsubara, H. Tamura, M. Takatsu, D. Yoo, B. Vatankhahghadim, H. Yamasaki, T. Miyazawa, S. Tsukamoto, Y. Watanabe, K. Takemoto, and A. Sheikholeslami, Complex, Intelligent, and Software Intensive Systems, edited by Leonard Barolli and Olivier Terzo (Springer International Publishing, Cham, 2018) 432–438.
  * [18] S. Tsukamoto, M. Takatsu, S. Matsubara, and H. Tamura, Fujitsu Sci. Tech. J 53, 8 (2017).
  * [19] S. Matsubara, M. Takatsu, T. Miyazawa, T. Shibasaki, Y. Watanabe, K. Takemoto, and H. Tamura, 2020 25th Asia and South Pacific Design Automation Conference (ASPDAC), pp. 667–672 (2020).
  * [20] H. Goto, K. Tatsumura, and A. R. Dixon, Science Advances 5 (2019).
  * [21] N. Mohseni, P. L. McMahon, and T. Byrnes, Nature Reviews Physics 4, 363, (2022).
  * [22] S. A. Cook, In Proceedings of the third annual ACM symposium on Theory of computing (STOC ’71), 151–158 (1971).
  * [23] L. Levin, Problems of Information Transmission 9, 115 (1973).
  * [24] M. Kowalsky, T. Albash, I. Hen, and D. A Lidar, Quantum Sci. Technol. 7, 025008 (2022).
  * [25] I. P. Gent and T. Walsh, ECAI 94, 105 (1994).
  * [26] F. Ivančić, Z. Yang, M. K. Ganai, et. al., Computer Aided Verification: 17th International Conference, CAV 2005, Edinburgh, Scotland, UK, July 6-10, 2005. Proceedings 17. Springer Berlin Heidelberg 2005.
  * [27] F. Ivančić, Z. Yang, M. Ganai, A. Gupta, P. Ashar, Theoretical Computer Science 404, 256 (2008).
  * [28] T. Grimm, D. Lettnin and M. Hübner Electronics, 7, 81 (2018).
  * [29] W. Gong, X. Zhou A survey of SAT solver AIP Conference Proceedings (Vol. 1836, No. 1). AIP Publishing
  * [30] A. Biere, M. Heule, H. van Maaren, T. Walsh:, Handbook of Satisfiability - Second Edition (Frontiers in Artificial Intelligence and Applications 336, IOS Press, 2021).
  * [31] V. Choi, arXiv:1004.2226 (2010).
  * [32] N. Chancellor, S. Zohren, P. A. Warburton, S. C. Benjamin, and S. Roberts, Sci. Rep. 6, 37107 (2016).
  * [33] J. Nüßlein, S. Zielinski, T. Gabor, C. Linnhoff-Popien, and S. Feld, Computational Science – ICCS 2023. ICCS 2023. Lecture Notes in Computer Science, vol 10477. Springer, Cham.
  * [34] S. Zielinski, J. Nüßlein, J. Stein, T. Gabor, C. Linnhoff-Popien, S. Feld, GECCO ’23 Companion: Proceedings of the Companion Conference on Genetic and Evolutionary Computation July 2023, Pages 2263–2271
  * [35] S. Zielinski, J. Nüßlein, J. Stein, T. Gabor, C. Linnhoff-Popien, S. Feld, Electronics 2023 12, 3492 (2023)
  * [36] J. D. Biamonte, Phys. Rev. A. 77, 052331 (2008).
  * [37] E. Boros, and A. Gruber, arXiv:1404.6538 (2014).
  * [38] M. Anthony, E. Boros, Y. Crama, and A. Gruber, Math. Program. 162, 115 (2017).
  * [39] N. Dattani, arXiv:1901.04405 (2019).
  * [40] A. Verma, M. Lewis, and G. Kochenberger, arXiv:2107.11695 (2021).
  * [41] Y. Crama, S. Elloumi, A. Lambert, and E. Rodriguez-Heck, hal-03795395 (2022).
  * [42] Fujitsu Services GmbH, Germany, Digital Annealer Tutorial Notebooks, (2023). Available at https://www.fujitsu.com/de/themes/digitalannealer/get-started/get-started-en.html

##  Appendix A Automatic Temperature Determination

A well adjusted search process for the Digital Annealer typically has two
phases. In the first phase the temperature is sufficiently high to instantly
escape from local minima. In this phase the random walk has progress in every
step. In the second phase temperatures should be low enough to produce waiting
cycles before escaping from a local minimum. Both phases and the transition
point from the first to the second phase can be controlled by the start and
end temperature. To estimate the dependency of the flip probability at a local
minimum from the temperature let $X$ be a minimum state and $X^{i}$ be the
state for which the bit with index $i$ is flipped:

$\displaystyle x^{i}_{j}=\begin{cases}x_{j}&\text{for }i\neq j\\\
1-x_{j}&\text{for }i=j\,.\end{cases}$ (11)

Let $T$ be the temperature and $E(X)$ and $E(X^{i})$ the energies of the
states. Since $X$ is a local minimum we can calculate the flip probability as
follows:

$\displaystyle
p_{\rm{flip}}=1-\prod_{i=0}^{N-1}\left(1-e^{-\frac{E(X_{i})-E(X)}{T}}\right)\,.$
(12)

Let $\Delta E(X)$ be a random variable defined on the set of minima of $E$ as
the mean of energy differences for all bit flips

$\displaystyle\Delta
E(X)=\frac{1}{N}\sum_{i=0}^{N-1}\left(E(X_{i})-E(X)\right)\,.$ (13)

As rough estimation we next replace all energy differences in Eq. (13) by the
expectation value of $\Delta E(X)$

$\displaystyle p_{\rm{flip}}\approx
1-\prod_{i=0}^{N-1}\left(1-e^{-\frac{\mathbb{E}(\Delta E(X))}{T}}\right)$ (14)

which leads to an estimation for a suitable temperature

$\displaystyle T\approx-\frac{\mathbb{E}(\Delta
E(X))}{\ln(1-\sqrt[N]{1-p_{\rm{flip}}})}\,.$ (15)

_Start Temperature_
In the beginning, waiting cycles in local minima should be avoided. Therefore,
we define a high flip probability, e.g., $p_{\rm{start}}=0.99$. According to
Eq. (15) we get for the start temperature

$\displaystyle T_{\rm{start}}=-\frac{\mathbb{E}(\Delta
E(X))}{\ln(1-\sqrt[N]{1-p_{\rm{start}}})}.$ (16)

_Transition Temperature_
For the transition point between the first and the second phase waiting cycles
in local minima should become more likely. Therefore, we define a lower flip
probability, e.g., $p_{\rm{trans}}=0.5$. According to Eq. (15) we get for the
transition temperature

$\displaystyle T_{\rm{trans}}=-\frac{\mathbb{E}(\Delta
E(X))}{\ln(1-\sqrt[N]{1-p_{\rm{trans}}})}\,.$ (17)

_End Temperature_
Let $I$ be the total number of iterations. Let the transition point be defined
by a part of iterations $0<\nu\leq 1$. Let $d:=1-\rm{decay}$ be the factor
multiplied for exponential temperature model in every step to the previous
temperature. Then we have $T_{\rm{start}}d^{I}=T_{\rm{end}}$ and
$T_{start}d^{\nu I}=T_{trans}$. The first equations gives
$d^{I}=T_{start}^{-1}T_{end}$ and using this to replace $d^{I}$ in the second
equation, we get
$T_{\rm{start}}T_{\rm{start}}^{-\nu}T_{\rm{end}}^{\nu}=T_{\rm{trans}}$.
Therefore
$T_{\rm{end}}=T_{\rm{start}}^{1-\frac{1}{\nu}}T_{\rm{trans}}^{-\frac{1}{\nu}}$
and with Eq. (17) we get

$\displaystyle
T_{\rm{end}}=T_{\rm{start}}^{1-\frac{1}{\nu}}\left(-\frac{\mathbb{E}(\Delta
E(X))}{\ln(1-\sqrt[N]{1-p_{\rm{trans}}})}\right)^{-\frac{1}{\nu}}\,.$ (18)

_Offset Parameter_
Finally we have to determine the dynamic energy offset that should help to
escape faster from local minima. To climb safely upwards from a local minimum,
we add only a fraction of the expected depth of the minimum in every step. If
we want to reach the rim in, e.g., $k:=10$ steps, then we define

$\displaystyle E_{\rm{dyn\\_off}}:=\frac{\mathbb{E}(\Delta E(X))}{k}\,.$ (19)