# The Perturbed Ferromagnetic Chain: A Tuneable Test of Hardness in the
Transverse-Field Ising Model

D.T. O’Connor  Department of Electronic and Electrical Engineering, University
College London, Gower Street, London, UK  London Centre for Nanotechnology,
Gordon Street, London, UK  L. Fry-Bouriaux  London Centre for Nanotechnology,
Gordon Street, London, UK  P.A. Warburton  Department of Electronic and
Electrical Engineering, University College London, Gower Street, London, UK
London Centre for Nanotechnology, Gordon Street, London, UK

(August 29, 2024)

###### Abstract

Quantum annealing in the transverse-field Ising model (TFIM) with open-system
dynamics is known to use thermally-assisted tunneling to drive computation.
However, it is still subject to debate whether quantum systems in the presence
of decoherence are more useful than those using classical dynamics to drive
computation. We contribute to this debate by introducing the perturbed
ferromagnetic chain (PFC), a chain of frustrated sub-systems where the degree
of frustration scales inversely with the perturbation introduced by a tunable
parameter. This gives us an easily embeddable gadget whereby problem hardness
can be tuned for systems of constant size. We outline the properties of the
PFC and compare classical spin-vector Monte Carlo (SVMC) variants with the
adiabatic quantum master equation. We demonstrate that SVMC methods get
trapped in the exponentially large first-excited-state manifold when solving
this frustrated problem, whereas evolution using quantum dynamics remains in
the lowest energy eigenstates. This results in significant differences in
ground state probability when using either classical or quantum annealing
dynamics in the TFIM.

##  I Introduction

Experimental validation of quantum processes and/or computational scaling
advantages in adiabatic quantum computation typically relies on the use of
artificial gadgets [1, 2, 3, 4, 5, 6], as these provide a way to demonstrate
quantum dynamics that can be exploited to enhance computation. A possible
framework to implement these gadgets is through the transverse-field Ising
model (TFIM), where the quantum dynamics are introduced through the addition
of a local non-commuting Hamiltonian to the system. If one interpolates
between the non-commuting Hamiltonians such that the system ends in the
computational basis of a problem space, then this is known as quantum
annealing [7, 8]. Such Hamiltonians that use this method typically have the
form

$\hat{H}(s)=-A(s)\sum_{i=1}^{N}\hat{\sigma}_{i}^{x}+B(s)\left[\sum_{i=1}^{N}h_{i}\hat{\sigma}_{i}^{z}+\sum_{i,j}J_{ij}\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}\right]\,,$
(1)

where $\hat{\sigma}^{x}$ and $\hat{\sigma}^{z}$ are the Pauli X and Z matrices
respectively, and the $n$-qubit problem is encoded in the biases, $h_{i}$, and
couplers, $J_{ij}$. The coefficients $A(s)$ and $B(s)$ are positive
monotonically decreasing and increasing functions in normalized time, $s$,
respectively. Typically the system is interpolated from the first term, which
is the transverse field component that has an easy-to-find ground state
solution at $s=0$, to the second term that encodes the problem of interest at
$s=1$.

This Hamiltonian has been extensively studied on experimental quantum
annealers for the past two decades, both within the context of combinatorial
optimization [9, 10, 11, 12, 13, 14, 15] and quantum simulation [16, 6, 17,
18, 19]. However, the quantumness of the dynamics used for computation on
quantum annealers is still subject to debate due to the prevalent quantum and
classical noise sources that can obscure coherent quantum processes [20, 21,
22, 23].

One such process is that of thermalization, which has been seen to aid
computation on quantum annealers where the anneal times are orders of
magnitude larger than the single-qubit decoherence time [24]. Improvements to
ground state probability can also be realised if we pause mid-anneal and allow
the system to thermalize near the minimum gap [25, 26, 27]. However, to what
extent thermalization is beneficial is still an active area of research, with
lower noise quantum annealers (i.e. those with lower thermalization rates)
being seen to improve tunneling ranges in local search (reverse) quantum
anneals [28]. It was also shown in Ref. [26] that the thermalization signature
of a quantum anneal can be replicated with spin-vector Monte Carlo (SVMC)
[29], which is a classical heuristic used to mimic the behaviour of physical
quantum annealers [30, 4]. Therefore it is important to ask where
thermalization with quantum dynamics can be useful computationally when
annealing in the TFIM.

In order to help answer this, we introduce the perturbed ferromagnetic chain
(PFC), a scaleable and tunable gadget that we use to differentiate between
SVMC and a system annealed under quantum dynamics. The PFC possesses the
qualities of a having a false minimum during an anneal (if the magnitude of
perturbation is small enough) and has an exponentially large degenerate first
excited state manifold in the computational basis. Additionally we can tune
the minimum gap energy, $\Delta_{10}$, with the perturbative parameter such
that it can be tuned through the value of the environmental temperature.

In the thermal regime, it is possible to compare the computational use of
thermally assisted quantum tunneling with classical mechanisms. The extent of
the former is analytically explored using the quantum adiabatic master
equation (AME) [31, 21], whereby we simulate quantum annealing using open
system dynamics. This provides a model of tunneling in a system that
experiences decoherence at a finite temperature, such that we can observe its
computational role when annealing with the PFC.

The format of the paper is as follows. We outline the definition of the PFC in
section II and give an overview of its properties, both classically at thermal
equilibrium and in the TFIM by inspecting the behaviour of both the quantum
and semi-classically approximated states. This will set the foundation as to
why this is a problem of interest when testing the quantumness of certain
methods in the TFIM. In section III we introduce the methods for simulating
both the quantum system using the AME and the classical system using SVMC. In
order to fully explore the dynamics of SVMC, we look at both SVMC and SVMC-TF
(SVMC with transverse-field dependent updates [26]) as well as introducing an
additional degree of freedom into both variants, which will allow for full
exploration of the Bloch sphere. Finally in section IV we present the
dynamical simulation results and demonstrate the effect that the false minimum
and the exponentially large first excited state manifold has on ground state
probability for both the AME and SVMC. We look at these effects with respect
to both the magnitude of perturbation and the PFC system size.

##  II Perturbed Ferromagnetic Chain

###  II.1 Classical Model

The PFC (Fig. 1) is a ferromagnetically coupled chain of frustrated sub-
systems each composed of two qubits. This system is similar to the cyclic
spin-gadgets used in Refs. [32, 33], but instead the gadget is made to be
acyclic and a perturbative offset, $d$, is introduced to break the degeneracy
of the ground state (at $d=0$) into a single ground state and
$2^{M}$-degenerate first excited state. The degree of frustration in the PFC
scales inversely with $d$. The PFC Hamiltonian is given by

$\displaystyle\frac{1}{R}\hat{H}_{P}=\hat{H}_{ss}-\sum^{M-1}_{i=1}\hat{\sigma}^{z}_{b,i}\hat{\sigma}^{z}_{b,i+1}\,,$
(2)
$\displaystyle\hat{H}_{ss}=\sum^{M}_{i=1}(1-d)\hat{\sigma}^{z}_{b,i}-\hat{\sigma}^{z}_{a,i}-\hat{\sigma}^{z}_{a,i}\hat{\sigma}^{z}_{b,i}\,,$

where $M$ is the number of subsystems, $R$ scales the energy of the problem,
and the magnitude of the perturbation is characterized by the parameter $d$.
The auxiliary and backbone qubits, depicted by the yellow and turquoise
circles shown in Fig. 1 respectively, are denoted by Pauli Z matrices
$\hat{\sigma}^{z}_{a}$ and $\hat{\sigma}^{z}_{b}$ respectively. The ground
state of this Hamiltonian is the $|0_{b}^{\otimes M},0_{a}^{\otimes M}\rangle$
(all up) state, and given that $d<(M-1)^{-1}$ the first excited state is a
degenerate manifold whose size grows exponentially in $M$. This manifold
always has the backbone qubits in the $|1_{b}\rangle^{\otimes M}$ (all down)
configuration, and the auxiliary qubits are iso-energetic with respect to
their spin state, creating a $2^{M}$-degenerate manifold of "floppy" auxiliary
qubits. This creates an energy gap between the ground state and exponential
manifold of $\Delta_{10}=2RMd$, where the states in the manifold have a
Hamming distance between $M$ and $2M$ from the ground state.

Figure 1: Illustration of the PFC Hamiltonian in general form. In yellow are
the auxiliary qubits with biases of $-R$, and in turquoise are the backbone
qubits with biases $R(1-d)$. The dark blue edges are the ferromagnetic
couplers of strength $-R$. The total number of qubits, $N$, in this problem is
$2M$, where $M$ denotes the number of 2-qubit sub-systems (indexed by $i$) in
the Hamiltonian. The system energy is scaled by $R$, and $d$ describes the
magnitude of perturbation. The properties of this model hold generally for
$M\geq 2$, $R>0$ and $1>d>0$.

Classically the PFC is exactly solvable via the transfer matrix method [34],
where at an inverse temperature, $\beta$, the partition function,
$\mathcal{Z}$, can be found in polynomial time. The partition function in
transfer matrix form is expressed as

$\displaystyle\mathcal{Z}=\mathbf{v}\mathbf{W}^{M-1}\mathbf{v}^{T}\,,$ (3)

where,

$\displaystyle\mathbf{v}=\begin{pmatrix}e^{\frac{1}{2}\beta
R\left(d+1\right)}&e^{\frac{1}{2}\beta R\left(d-3\right)}&e^{\frac{1}{2}\beta
R\left(1-d\right)}&e^{\frac{1}{2}\beta R\left(1-d\right)}\\\ \end{pmatrix}$
(4)

handles the boundary subsystems of the chain, and

$\displaystyle\mathbf{W}=\begin{bmatrix}e^{\beta R\left(d+2\right)}&e^{\beta
Rd}&1&1\\\ e^{\beta Rd}&e^{\beta R\left(d-2\right)}&e^{-2\beta R}&e^{-2\beta
R}\\\ 1&e^{-2\beta R}&e^{\beta R\left(2-d\right)}&e^{\beta
R\left(2-d\right)}\\\ 1&e^{-2\beta R}&e^{\beta R\left(2-d\right)}&e^{\beta
R\left(2-d\right)}\end{bmatrix},$ (5)

handles the inner subsystems of the chain. We can then find the magnetization
of a subsystem at thermal equilibrium using

$\displaystyle\langle\sigma^{z}_{i}\rangle=\frac{1}{\mathcal{Z}}\left[\mathbf{v}\mathbf{W}^{i-1}\sigma^{z}_{i}\mathbf{W}^{M-i}\mathbf{v}^{T}\right]\,,$
(6)

where
$\sigma^{z}_{i}=\frac{1}{2}\left(\hat{\sigma}^{z}_{a,i}+\hat{\sigma}^{z}_{b,i}\right)$.
The average magnetization of the PFC is then an average over all
contributions,
$\langle\sigma^{z}\rangle=\frac{1}{M}\sum_{i}^{M}\langle\sigma^{z}_{i}\rangle$.

Additionally, the free-energy of the PFC can be derived from the transfer
matrix, which can then be used to derive further thermodynamic properties. The
free-energy, $F$, is defined as

$F=-\lim_{M\rightarrow\infty}\frac{1}{\beta M}\ln\mathcal{Z}\,.$ (7)

We show in Appendix A that in the limit of $M\rightarrow\infty$, that
$\mathcal{Z}^{1/M}\rightarrow\lambda_{1}$, where $\lambda_{1}$ is the spectral
radius (largest absolute eigenvalue) of $\mathbf{W}$. The free-energy of the
PFC is found to be

$\begin{split}F&=-\frac{1}{\beta}\ln\Bigg{(}e^{2\beta R}\cosh{\beta
Rd}+\cosh{\beta R(2-d)}+\\\ &\sqrt{\left(e^{2\beta R}\cosh{\beta
Rd}+\cosh{\beta R(2-d)}\right)^{2}-4\sinh{4\beta R}}\Bigg{)}\,,\end{split}$
(8)

which is always real, continuous and finite given that the original parameter
constraints of the PFC are met.

###  II.2 Transverse-Field Ising Model

(a)

(b)

Figure 2: (a) Intensity plot of the average qubit magnetization in the
instantaneous ground state for a PFC in the presence of a transverse-field
(Eq. (9)) with $M=2$ and $R=1.0$. The red line shows the boundary between the
quantum paramagnetic and negative magnetization phases. The white dashed line
indicates the position of the minimum gap. (b) Cross-sections of (a) showing
the average magnetization during an anneal. Figure 3: Plots of the semi-
classical energy potential (top row), and the energy along the hyper-plane
(white dashed line) passing through the landscape (bottom row) at the
specified values of normalized time, $s$. This potential is for a PFC with
$M=2$, $R=1.0$ and $d=0.09$. The red marker indicates the global minimum of
the landscapes. The backbone and auxiliary spins in the PFC are parameterized
by angles $\theta_{b}$ and $\theta_{a}$ respectively.

Translating the PFC into the TFIM involves the addition of a non-commuting
transverse field term, composed of $\hat{\sigma}^{x}$ operators, which
introduces the driver of quantum fluctuations that can potentially be used to
aid computation [7, 8]. The TFIM Hamiltonian of the PFC is given by

$\hat{H}(s)=-A(s)\sum_{j=1}^{N}\hat{\sigma}_{j}^{x}+B(s)\hat{H}_{P}\,,$ (9)

where the classical PFC is encoded into $\hat{H}_{P}$ (Eq. (2)). The
coefficients are taken to be $A(s)=3(1-s)$ and $B(s)=3s$ throughout this work,
where $s$ is the normalised annealing time $s=t/t_{\textrm{anneal}}$.

For sufficiently small values of $d$, the ground state of this Hamiltonian is
seen to undergo a quantum phase transition, illustrated in Fig. 2 by the
change in average qubit magnetization from negative to positive phases. The
average qubit magnetization is defined as

$\langle\sigma^{z}\rangle=\frac{1}{N}\sum_{j=1}^{N}\langle
E_{0}(s)|\hat{\sigma}^{z}_{j}|E_{0}(s)\rangle\,,$ (10)

where $|E_{0}(s)\rangle$ is the instantaneous ground state from the
diagonalized Hamiltonian of Eq. (9) at some value of $s$, and $N$ is the
number of qubits in the PFC. The formation of the negative phase before the
minimum gap is indicative of the ground state qubits becoming magnetized to
resemble the exponentially large degenerate first-excited-state manifold
(further illustrated in sections II.3 and IV). After passing through the
minimum gap, the instantaneous ground state enters the positive phase and then
goes on to finish in the computational ground state. It must be noted that
even if the $d<(M-1)^{-1}$ condition is broken such that the exponentially
degenerate manifold is no longer the computational first-excited-state, for
small $d$ the instantaneous ground state maintains its resemblance to the
exponentially degenerate manifold before the minimum gap.

The tuning of $d$ scales the size of the gap for all $M$ in the TFIM, and
therefore can be used as a tunable hardness parameter for the PFC. Although
scaling with $M$ is numerically seen to exponentially scale the minimum gap
size and therefore hardness, it can lead to intractable computational times
for some of the simulation methods explored in the results, making $d$ a more
desirable tunable hardness parameter.

In summary, the PFC becomes hard in the TFIM for small values of $d$ due to
the presence of a quantum phase transition, and for large values of $M$ where
the minimum gap exponentially reduces in size. In the next section we perform
a semi-classical analysis to show that when the PFC is translated into to the
TFIM a false minimum exists. The interplay between the false and true minima
is particularly prominent when annealing through the region where the ground
state is in the negative magnetization phase shown in Fig. 2.

###  II.3 Semi-Classical Analysis

To further explore the behaviour of the PFC in the TFIM, a semi-classical
approximation can be made by using the spin-coherent Ansatz [35]

$\left|\theta\right\rangle=\bigotimes^{N}_{j=1}\cos\left(\frac{\theta_{j}}{2}\right)|0\rangle+e^{i\phi_{j}}\sin\left(\frac{\theta_{j}}{2}\right)|1\rangle$
(11)

to calculate the semi-classical effective potential landscape as a function of
$s$. The magnetization expectation values of the auxiliary qubits are almost
identical, and the same is true for the backbone qubits. We can therefore
approximate the states of the PFC in the spin-coherent Ansatz as

$\displaystyle\left|\theta_{a},\theta_{b}\right\rangle$
$\displaystyle=\begin{aligned}
&\left[\bigotimes^{M}_{i=1}\cos\left(\frac{\theta_{a}}{2}\right)|0\rangle+\sin\left(\frac{\theta_{a}}{2}\right)|1\rangle\right]\otimes&\\\
&\left[\bigotimes^{M}_{i=1}\cos\left(\frac{\theta_{b}}{2}\right)|0\rangle+\sin\left(\frac{\theta_{b}}{2}\right)|1\rangle\right]\,.\end{aligned}$
(12)

Here, $\theta_{a}$ and $\theta_{b}$ are the angles of the states in the XZ-
plane of the Bloch spheres for all of the auxiliary and backbone qubits
respectively. Here we assume that the azimuthal angle, $\phi_{j}$ is equal to
zero. The semi-classical potential is then given by

$V_{SC}(s,\theta_{a},\theta_{b})=\left\langle\theta_{a},\theta_{b}\right|\hat{H}(s)\left|\theta_{a},\theta_{b}\right\rangle.$
(13)

The visual representation of the potential at various stages of an anneal
(Fig. 3) shows the PFC initially taking a path to the first excited states
($\theta_{a}\in[-\pi,\pi],\theta_{b}=\pi$). This is then followed by a
discontinuous change in the position of the global energy minimum about the
minimum gap (at $s=0.841$) to the computational ground state
($\theta_{a}=0,\theta_{b}=0$). By taking a hyper-plane that passes through the
global minimum (and the local minimum where applicable), it is clear that as
we evolve from a uni-modal to a bi-modal potential the all down
($\theta_{b}=\pi$) configuration of the backbone qubits is energetically
preferable until the minimum gap is traversed.

Additionally, the computational ground state is energetically isolated from
the low-energy excited states, meaning that further dynamical evolution is
still needed to reach the ground state after the minimum gap. If the system is
evolved under an adiabatic, coherent regime [36] then the dynamical process is
quantum tunneling. In a classical model (like SVMC) we can only use thermal
excitations to traverse these energy barriers.

If tunneling were to occur in the instantaneous ground state, then this would
result in delocalization about the bistable potential. Using the trace-norm
distance,

$D(s,\theta_{a},\theta_{b})=\sqrt{1-\left|\langle
E_{0}(s)|\theta_{a},\theta_{b}\rangle\right|^{2}}\,,$ (14)

we can quantify the distance between the instantaneous ground state and the
spin-coherent Ansatz, to determine the extent to which the Ansatz accurately
describes the instantaneous ground state. In Fig. 4 we show the trace-norm
distance for a four-qubit ($M=2$) instance of the PFC in the vicinity of the
minimum gap. This is the same instance whose potential landscape is shown in
Fig. 3 and whose minimum gap occurs at $s=0.841$. At $s=0.835$ the trace-norm
distance has a global minimum whose location in $(\theta_{a},\theta_{b})$
space closely corresponds to the global minimum of the semi-classical
potential (i.e. near $\theta_{b}=\pi$). Nevertheless there is a local minimum
in the trace-norm distance which extends along the indicated hyper-plane,
showing that, prior to the minimum gap, tunneling enables a finite probability
amplitude in the local minimum of the potential landscape near $\theta_{b}=0$.
After the minimum gap, as shown for the case $s=0.845$ in the right panel of
Fig. 4, the global minimum of the trace-norm distance now closely corresponds
to the global minimum of the semi-classical potential near $\theta_{b}=0$.
This minimum continuously evolves into the global minimum of the problem
Hamiltonian at $s=1$ as shown in Fig. 3.

Figure 4: The trace-norm distance between the spin-coherent state and the
instantaneous ground state (Eq. (14)) in the vicinity of the minimum gap at
$s=0.841$ for a PFC with $M=2$, $R=1.0$ and $d=0.09$. The backbone and
auxiliary spins in the PFC are characterised by angles $\theta_{b}$ and
$\theta_{a}$ respectively. The red cross marker indicates the global minimum
of the trace-norm landscapes and the red circle indicates the minimum of the
potential landscape in Fig. 3. The hyper-plane (white dashed line) from Fig. 3
is also plotted.

Using semi-classical analysis we have shown the existence of a false global
minimum before the minimum gap, when the value of $d$ is small enough. Under
quantum evolution, the transition from the false minimum to the true minimum
exploits tunnelling, and this is visualized by measuring the trace-norm
distance between spin-coherent states and the instantaneous ground state to
show delocalization across the potential barrier (Fig. 4). Under classical
evolution, thermal excitation of multiple qubits is needed to traverse the
barrier to reach the computational ground state. Additionally the manifold
along $\theta_{b}=\pi$ (corresponding to the exponentially large computational
first excited state manifold) is equally accessible under classical dynamics.
Therefore classical algorithms that explore this energy landscape, such as
SVMC, can remain in this manifold instead of reaching the ground state. We
explore in section IV the extent to which this hinders computation in both
classical and quantum evolution.

##  III Methods

###  III.1 Spin-Vector Monte Carlo

The spin-vector Monte Carlo algorithm [29] is an algorithm that replaces the
Pauli matrices in the Hamiltonian of Eq. (9) with $O(2)$ rotors in the XZ-
plane of the Bloch-sphere. The energy function we are to minimize using
Metropolis-Hastings updates becomes

$\begin{split}\mathcal{E}(s)&=-A(s)\sum^{N}_{j=1}\sin\theta_{j}+\\\
&B(s)\left[\sum^{N}_{j=1}h_{j}\cos\theta_{j}+\sum_{\langle
j,k\rangle}J_{jk}\cos\theta_{j}\cos\theta_{k}\right]\,,\end{split}$ (15)

where $A(s)$ and $B(s)$ are the same schedule functions introduced in Eq. (1).
The SVMC algorithm attempts to update all rotor angles, $\theta$, in every
sweep (i.e. every increment of time). The method of update can be described in
two ways:

$\displaystyle\theta_{j}^{t}\in[0,\pi],\quad\theta_{j}^{t+1}\sim\textrm{Uniform}(0,\pi)\,,$
(16)
$\displaystyle\theta_{j}^{t}\in[0,\pi],\quad\theta_{j}^{t+1}=\theta_{j}^{t}+\min\left(\frac{A(s)}{B(s)},1\right)u,$
$\displaystyle\quad u\sim\textrm{Uniform}(-\pi,\pi)\,.$ (17)

The traditional method of update in SVMC is described by Eq. (16), whereby the
new angle is a sample from a uniform distribution from zero to $\pi$. A more
recent update method used to capture additional annealing artifacts such as
freeze-out is shown in Eq. (17), which we refer to as SVMC-TF [26]. In this
latter version, the freedom of the rotor movement in an update is proportional
to the relative magnitude of the transverse field that drives the dynamics.

We also consider another variant of SVMC and SVMC-TF, whereby the dynamical
restriction of only operating in the XZ-plane is removed by including the
azimuthal angle $\phi_{j}$ in Eq. (12), such that SVMC now has access to the
entire Bloch sphere. We will refer to these variants as spherical-SVMC and
spherical-SVMC-TF. This coordinate extension does not affect the Z (polar)
components of the energy function. It does however affect the transverse field
component, such that the new energy function becomes

$\begin{split}\mathcal{E}(s)&=-A(s)\sum_{j=1}^{N}\cos\phi_{j}\sin\theta_{j}+\\\
&B(s)\left[\sum_{j=1}^{N}h_{j}\cos\theta_{j}+\sum_{\langle
j,k\rangle}J_{jk}\cos\theta_{j}\cos\theta_{k}\right]\,,\end{split}$ (18)

where the azimuthal angle, $\phi_{j}\in[-\pi,\pi]$, is also updated in the
same way as the polar angle, $\theta_{j}$.

Throughout the rest of this paper, we look at both the simplest case, SVMC,
and the more complex spherical-SVMC-TF for the comparative experiments, with
the other variants included in appendix B for completeness. For all variants,
we also take the annealing functions to be $A(s)=3(1-s)$ and $B(s)=3s$, where
we start the annealing from $s=0$ to $s=1$ at a temperature of $12$ mK. All
algorithms update spins individually in a randomly permuted order, thus cannot
capture any simultaneous multi-qubit moves unlike those that may occur in a
system evolved using quantum dynamics.

###  III.2 Quantum Dynamics

To simulate the quantum evolution of the PFC, we look at the dynamics in both
closed and open systems using the von-Neumann equation and adiabatic master
equation [31, 21] (AME) respectively. For the closed system simulations, the
von-Neumann equation takes the form,

$\hbar\frac{\partial}{\partial
t}\hat{\rho}(t)=-i\left[\hat{H}(t),\hat{\rho}(t)\right]\,,$ (19)

where $\hat{\rho}(t)$ is the density matrix, and $t$ is related to normalized
time, $s$, by $s=t/t_{\textrm{anneal}}$. The initial state is the pure ground
state of the system at $s=0$, which when using Eq. (9) is
$\hat{\rho}(0)=|+\rangle\langle+|$. To solve the von-Neumann equation and the
AME we have used the Hamiltonian Open Quantum System Toolkit [37] (HOQST), a
simulation library written for the Julia language.

The form of the AME we have chosen for the open system simulations is one that
uses parameters similar to those used in theoretical studies of the
experimental annealing hardware [33, 38]. We will not define the AME in its
entirety here (see Refs. [37, 31]), but we use a Davies-style AME that is
valid in the weak-coupling limit. It is also necessary to make assumptions on
the model of decoherence used, where we assume that all qubits are coupled
equally to the bath independently. All qubits experience decoherence by
dephasing, and the bath spectral density takes the form of a Bosonic Ohmic
bath,

$\gamma(\omega)=2\pi\eta
g^{2}\frac{\omega\exp\left(-|\omega|/\omega_{c}\right)}{1-\exp\left(-\beta\omega\right)}\,,$
(20)

where $\beta=(k_{B}T)^{-1}$ is the inverse temperature, $\omega_{c}$ is the
cutoff frequency, $\eta g^{2}$ is the dimensionless bath coupling strength and
$k_{B}$ is the Boltzmann constant. Throughout the rest of the paper we specify
the bath parameters to be $T=12$ mK, $\omega_{c}=4$ GHz and $\eta
g^{2}=10^{-3}$.

The decoherence by dephasing then manifests itself through the time-dependent
Lindblad operators

$\hat{L}_{j,\omega_{kl}}(t)=\langle
E_{l}(t)|\sigma^{z}_{j}|E_{k}(t)\rangle|E_{l}(t)\rangle\langle E_{k}(t)|\,.$
(21)

This describes how the $j^{th}$ qubit couples to environment according to
energy gaps, $\omega_{kl}=E_{k}-E_{l}$, between the instantaneous energy
eigenstates, $|E_{k}(t)\rangle$ and $|E_{l}(t)\rangle$ of the system
Hamiltonian (Eq. (9)).

We use this form of the AME since it models thermally assisted adiabatic
quantum computation near the adiabatic limit. Assuming that most of the ground
state population is lost to the first excited state after passing through the
minimum gap, the re-population of the ground state via thermal relaxation can
be related to the transition rate via

$\Gamma_{1\rightarrow 0}(t)\propto\gamma_{1\rightarrow
0}(t)=\gamma(\omega_{10}(t))\sum_{j}\left|\langle
E_{0}(t)|\sigma^{z}_{j}|E_{1}(t)\rangle\right|^{2}\,.$ (22)

Here the temperature dependence of the transition rate is introduced by the
Ohmic spectral density function, $\gamma$ in Eq. (20).

However, when the gap, $\omega_{10}$, is sufficiently small, the weak-coupling
assumption in the AME starts to break down. This is the case for hard PFC
instances since the gap can be very small relative to the bath temperature. In
such a regime the energy levels become broadened due to the stronger coupling
to the bath, such that the discrete energy levels should emulate a more
continuous potential, similar to the semi-classical picture. Therefore,
despite the AME not being able to describe these strong-coupling regimes as
accurately as more sophisticated models like the Redfield equation, it serves
as a reasonable approximation of an open-system model of the PFC.

##  IV Dynamical Simulations

Figure 5: Probability of being in the ground state for both SVMC (solid lines)
and spherical-SVMC-TF (dotted lines) as the system scales in size, $M$, for a
PFC with $d=0.1$. Here, the SVMC and spherical-SVMC-TF probabilities are found
from 20,000 samples, which we repeat 50 times and bootstrap to find the median
and $95\%$ confidence intervals for the data point and error bars
respectively. Figure 6: State Probability at $s=0.83$
($=t/t_{\textrm{anneal}}$), just after passing through the minimum gap
($s=0.8227$) for a) the AME with $t_{\textrm{anneal}}=200ns$ and b) spherical-
SVMC-TF with $t_{\textrm{anneal}}=10,000$ sweeps. We measure a PFC with $M=3$
and $d=0.1$. For the AME we measure the probability of being in the $j^{th}$
instantaneous state, $E_{j}$, and for spherical-SVMC-TF we take a classical
measurement of being in either the ground state or any of the first excited
states. The spherical-SVMC-TF probabilities are found from 20,000 samples, of
which we repeat 50 times and bootstrap to find a median and $95\%$ confidence
intervals for the data point and error bars respectively. Figure 7: Plots of
state probability for being in either the ground state (solid lines) or any of
the $2^{M}$-degenerate $1^{st}$ excited states (dashed lines) at the end of an
anneal. A PFC with system size $M=3$ was evolved using quantum (top row) and
classical (bottom row) dynamics. The black dashed line indicates the
probability with random guessing i.e. $1/64$. The quantum simulations are
plotted against anneal time in $ns$, whilst the SVMC simulations are plotted
against the number of sweeps. The closed and open system dynamics are evolved
according to the von-Neumann and adiabatic master equation respectively (see
section III.2). The SVMC and spherical-SVMC-TF probabilities are found from
20,000 samples, of which we repeat 50 times and bootstrap to find a median and
$95\%$ confidence intervals for the data point and error bars respectively.

We begin the dynamical analysis of the PFC by observing the performance of the
SVMC variants when scaling in $M$. The systems chosen meet the $d<(M-1)^{-1}$
condition such that the first excited state is the exponentially large
manifold which is at least a Hamming distance of $M$ away from the ground
state. The combined effects of an exponentially scaling gap and manifold are
observed in Fig. 5 by measuring the ground state probability at the end of the
anneal with respect to the number of incremental sweeps used in both SVMC and
spherical-SVMC-TF. Typically we expect an increasing number of sweeps to
correspond to an increasing ground state probability, but here there are three
distinct regimes when annealing the PFC. For low sweep numbers, where the
semi-classical potential is evolved in large steps, we have a relatively high
ground state probability as the false minimum is not well resolved but SVMC
still guides the spin-vector to the low-energy states.

For medium sweep numbers we see reduction in ground state probability, caused
by the quasi-continuous evolution of the semi-classical potential now leading
the SVMC algorithm to the false minimum. This guides SVMC to the
$\theta_{b}=\pi$ manifold (Fig. 3) corresponding to the degenerate
computational first-excited-states, causing SVMC to spread out into this
manifold and into states that are potentially further in Hamming distance from
the computational ground state (see Fig. 6 for further evidence of this).
Finally, for high sweep numbers SVMC starts to thermally equilibrate and
ground state probability begins to return.

To confirm the detrimental role that the exponential manifold has on SVMC, the
state probabilities just after the minimum gap are measured in Fig. 6 for both
spherical-SVMC-TF (for $10,000$ sweeps) and a system evolved using AME (for
$200ns$). The probability distribution was measured at $s=0.83$ for a PFC
($M=3$, $R=1.0$, $d=0.1$) with a minimum gap at $s=0.8227$. It can be seen
that the probability density spreads into the computational first excited
state manifold when evolving using spherical-SVMC-TF, which is caused by
following the false minimum and accessing the $\theta_{b}=\pi$ manifold using
classical dynamics. However, a system evolved using the AME remains in the
lowest instantaneous eigenstates after the minimum gap, where $E_{0}$ also has
a large overlap with the computational ground state, $\langle
E_{0}(s=0.83)|0^{\otimes N}\rangle=0.98$. After this point in the anneal, both
the AME and spherical-SVMC-TF experience freeze-out, which prevents any more
dynamical evolution that could affect ground state probability. This is
illustrated by the ground state probability of the AME (spherical-SVMC-TF)
evolution at $s=0.83$ being $\sim 0.50$ $(\sim 3.4\times 10^{-3})$, compared
to that at $s=1$ being $\sim 0.51$ $(\sim 3.4\times 10^{-3})$ \- see Fig. 7.
It is also worth noting that we measure in the computational basis for SVMC
since it is a classical algorithm with no other analogous discrete states to
compare against the instantaneous states used by the AME.

Finally, we explore how the tunable hardness parameter, $d$, affects the PFC
in both quantum and SVMC simulations for $M=3$, in Fig. 7. We measure the
probability at the end of the anneal of being in the ground state as well as
any of the first excited states for different annealing durations. The value
of $d$ also determines the size of the minimum gap, such that we span $d$ to
capture various regimes at a fixed system temperature of $12$ mK. At $d\sim
0.227$ we have a minimum gap that approximately equals the system temperature.

The form of SVMC and spherical-SVMC-TF when scaling in $d$ in Fig. 7 is
similar to what is also seen in Fig. 5 when scaling in $M$. For all sweeps,
the PFC becomes harder as $d$ becomes smaller, and SVMC and spherical-SVMC-TF
preferably anneal to the first excited state manifold (following the false
minimum). Additionally, spherical-SVMC-TF consistently out-performs SVMC for
our hardest problems ($d\leq 0.15$), something which is further discussed in
appendix B.

The probabilities at the end of a closed-system quantum anneal (for time
$t_{\textrm{anneal}}$) are of a similar form, with a bump for short anneal
times (the diabatic bump [39]) and then reaching the adiabatic limit (where
the ground state probability tends towards 1) at longer anneal times (e.g.
$t_{\textrm{anneal}}\simeq 20ns$ for $d=0.3$). The closed system dynamics
describe quantum evolution at a temperature of $0K$ with no dephasing, and
therefore requires a run time of $t_{\textrm{anneal}}=O(1/\Delta_{10}^{2})$ to
run adiabatically. For example, at $d=0.05$ the approximate adiabatic run time
is $\sim 163\mu s$ and we therefore see a probability $\sim 0$ due to the
short run times (see appendix B for a log-scale plot of the closed-system
quantum anneal).

However, for the open system simulations we see a non-zero ground state
probability at $d=0.05$ for $t_{\textrm{anneal}}=200ns$. This can be
attributed to thermally assisted quantum transitions, whereby the relaxation
rate (Eq. (22)) is non-zero about the minimum gap (Fig. 8) and can return
probability to the ground state from the first excited state. This can occur
because the system is sub-thermal [27] (i.e. the ground-state probability is
less than that at thermal equilibrium) immediately after the minimum gap, and
the energy gap is still small enough to allow significant thermalization from
the first excited state to the ground state. Additionally, this transition
involves all backbone qubits changing their magnetization simultaneously.

Figure 8: Plot of the evolution of ground state probability from the AME
(black solid) and the Gibbs state (black dotted), as well as the transition
rate (blue) (Eq. 22). A PFC of $M=3$, $R=1.0$, and $d=0.05$ was simulated at a
system temperature of 12 mK. The AME was evolved for
$t_{\textrm{anneal}}=200ns$ and the minimum gap occurs at $s=0.9059$.

Therefore thermalization in a quantum system is seen to be of some
computational use, as has been seen in other literature [24, 25, 27, 26, 28].
However, we see that on the time-scales tested that thermalization in
classical SVMC plays a far less significant role, and results in a ground
state probability (after $10,000$ sweeps) two orders of magnitude lower than a
$100ns$ AME evolved anneal. The marked difference between the two types of
dynamical simulations highlight the effect of the minimum gap and the
exponential manifold on ground state probability, making it a gadget of
interest for when differentiating between quantum and classical evolutions in
the TFIM. For a better contrast between dynamics, annealing larger versions of
the PFC ($M>10$) would result in negligible statistical contribution from
random state sampling, and an extremely large first excited state manifold
that would likely result in SVMC failing to find the ground state. Open system
simulations of this kind are unfeasible, and would most likely only be
realised on experimental quantum hardware.

##  V Conclusion

In this work, we have introduced the perturbed ferromagnetic chain (PFC), a
gadget with an exponentially large first-excited-state manifold and an
isolated ground state, whereby problem hardness and frustration is tuned by
the perturbative parameter, $d$. When annealed in the transverse field Ising
model (TFIM), the PFC develops computationally hard characteristics such as an
exponentially small minimum gap (in $N$), a quantum phase transition, and a
false minimum.

The evolution of the PFC in the TFIM was assessed with quantum dynamics using
the adiabatic master equation (AME), and classically using both spin-vector
Monte Carlo (SVMC) and spherical-SVMC-TF (see section II.3 for more
information). For quasi-continuous evolution of the PFC with the SVMC methods,
the false minimum is followed to the computational first-excited-state
manifold. This results in probable transitions to other low energy states
further in Hamming distance away from the computational ground state (Fig. 6),
and therefore reduces the probability of reaching the ground state
significantly. This is compounded by increasing problem size (Fig. 5) and by
tuning $d$ (Fig. 7).

For a PFC evolved using the AME, the system mostly remains in the lowest two
eigenstates (Fig. 6) instead of accessing the exponential manifold, such that
a $100ns$ open-quantum-system anneal results in a ground state probability
two-orders of magnitude larger than a $10,000$ sweep spherical-SVMC-TF anneal
for the hardest comparative problem simulated (Fig. 7). The AME evolution
permits thermalization to the ground state (Fig. 8) in time-scales too short
for adiabatic evolution, indicating that thermalization aids computation of
the PFC in the TFIM.

The distinct differences seen between classical and quantum evolutions
therefore makes the PFC a useful gadget in differentiating dynamical systems,
something which other gadgets cannot always exhibit [4, 30, 26]. Comparative
tests at larger PFC system sizes of $M>10$ would also allow for better
differentiation and insight into the computational extent of thermalization in
quantum and classical evolutions. However, realizing open quantum system
simulations of this scale is intractable, and would therefore require further
approximations or an experimental implementation with quantum hardware.

## Acknowledgements

We gratefully acknowledge Catherine White, Gioele Consani, Robert Banks and
Nedeen Alsharif for useful comments and insightful discussion. The research is
based upon work (partially) supported by the Office of the Director of
National Intelligence (ODNI), Intelligence Advanced Research Projects Activity
(IARPA), and the Defense Advanced Research Projects Agency (DARPA), via the
U.S. Army Research Office Contract No. W911NF-17-C-0050. The views and
conclusions contained herein are those of the authors and should not be
interpreted as necessarily representing the official policies or endorsements,
either expressed or implied, of the ODNI, IARPA, DARPA, or the U.S.
Government. The U.S. Government is authorized to reproduce and distribute
reprints for Governmental purposes notwithstanding any copyright annotation
thereon. D. T. O’C. thanks EPSRC and BT for additional support.

## References

  * Nagaj  _et al._ [2012] D. Nagaj, R. D. Somma, and M. Kieferova, Quantum Speedup by Quantum Annealing, Physical Review Letters 109, 050501 (2012), arXiv:1202.6257 .
  * Denchev  _et al._ [2016] V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Babbush, V. Smelyanskiy, J. Martinis, and H. Neven, What is the Computational Value of Finite-Range Tunneling?, Physical Review X 6, 031015 (2016).
  * Boixo  _et al._ [2016] S. Boixo, V. N. Smelyanskiy, A. Shabani, S. V. Isakov, M. Dykman, V. S. Denchev, M. H. Amin, A. Y. Smirnov, M. Mohseni, and H. Neven, Computational multiqubit tunnelling in programmable quantum annealers, Nature Communications 7, 10327 (2016).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Demonstration of a scaling advantage for a quantum annealer over simulated annealing, Physical Review X 8, 031016 (2018), arXiv:1705.07452 .
  * Mandrà and Katzgraber [2018] S. Mandrà and H. G. Katzgraber, A deceptive step towards quantum speedup detection, Quantum Science and Technology 3, 04LT01 (2018).
  * King  _et al._ [2021] A. D. King, J. Raymond, T. Lanting, S. V. Isakov, M. Mohseni, G. Poulin-Lamarre, S. Ejtemaee, W. Bernoudy, I. Ozfidan, A. Y. Smirnov, M. Reis, F. Altomare, M. Babcock, C. Baron, A. J. Berkley, K. Boothby, P. I. Bunyk, H. Christiani, C. Enderud, B. Evert, R. Harris, E. Hoskinson, S. Huang, K. Jooya, A. Khodabandelou, N. Ladizinsky, R. Li, P. A. Lott, A. J. R. MacDonald, D. Marsden, G. Marsden, T. Medina, R. Molavi, R. Neufeld, M. Norouzpour, T. Oh, I. Pavlov, I. Perminov, T. Prescott, C. Rich, Y. Sato, B. Sheldan, G. Sterling, L. J. Swenson, N. Tsai, M. H. Volkmann, J. D. Whittaker, W. Wilkinson, J. Yao, H. Neven, J. P. Hilton, E. Ladizinsky, M. W. Johnson, and M. H. Amin, Scaling advantage over path-integral Monte Carlo in quantum simulation of geometrically frustrated magnets, Nature Communications 12, 1113 (2021).
  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse Ising model, Physical Review E 58, 5355 (1998).
  * Farhi  _et al._ [2001] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, A Quantum Adiabatic Evolution Algorithm Applied to Random Instances of an NP-Complete Problem, Science 292, 472 (2001).
  * Venturelli and Kondratyev [2019] D. Venturelli and A. Kondratyev, Reverse Quantum Annealing Approach to Portfolio Optimization Problems, Quantum Machine Intelligence 1, 17 (2019), arXiv:1810.08584 .
  * Stollenwerk  _et al._ [2018] T. Stollenwerk, E. Lobe, and M. Jung, Flight Gate Assignment with a Quantum Annealer, arXiv:1811.09465 [quant-ph] (2018), arXiv:1811.09465 [quant-ph] .
  * Kim  _et al._ [2019] M. Kim, D. Venturelli, and K. Jamieson, Leveraging Quantum Annealing for Large MIMO Processing in Centralized Radio Access Networks, Proceedings of the ACM Special Interest Group on Data Communication , 241 (2019), arXiv:2001.04014 .
  * Inoue  _et al._ [2021] D. Inoue, A. Okada, T. Matsumori, K. Aihara, and H. Yoshida, Traffic signal optimization on a square lattice with quantum annealing, Scientific Reports 11, 3303 (2021).
  * Kitai  _et al._ [2020] K. Kitai, J. Guo, S. Ju, S. Tanaka, K. Tsuda, J. Shiomi, and R. Tamura, Designing metamaterials with quantum annealing and factorization machines, Physical Review Research 2, 013319 (2020).
  * Asproni  _et al._ [2020] L. Asproni, D. Caputo, B. Silva, G. Fazzi, and M. Magagnini, Accuracy and minor embedding in subqubo decomposition with fully connected large problems: A case study about the number partitioning problem, Quantum Machine Intelligence 2, 1 (2020).
  * Neukart  _et al._ [2017] F. Neukart, G. Compostella, C. Seidel, D. von Dollen, S. Yarkoni, and B. Parney, Traffic Flow Optimization Using a Quantum Annealer, Frontiers in ICT 4, 10.3389/fict.2017.00029 (2017).
  * King  _et al._ [2018] A. D. King, J. Carrasquilla, J. Raymond, I. Ozfidan, E. Andriyash, A. Berkley, M. Reis, T. Lanting, R. Harris, F. Altomare, K. Boothby, P. I. Bunyk, C. Enderud, A. Fréchette, E. Hoskinson, N. Ladizinsky, T. Oh, G. Poulin-Lamarre, C. Rich, Y. Sato, A. Y. Smirnov, L. J. Swenson, M. H. Volkmann, J. Whittaker, J. Yao, E. Ladizinsky, M. W. Johnson, J. Hilton, and M. H. Amin, Observation of topological phenomena in a programmable lattice of 1,800 qubits, Nature 560, 456 (2018).
  * Bando  _et al._ [2020] Y. Bando, Y. Susa, H. Oshiyama, N. Shibata, M. Ohzeki, F. J. Gómez-Ruiz, D. A. Lidar, S. Suzuki, A. del Campo, and H. Nishimori, Probing the universality of topological defect formation in a quantum annealer: Kibble-Zurek mechanism and beyond, Physical Review Research 2, 033369 (2020).
  * Harris  _et al._ [2008] R. Harris, M. W. Johnson, S. Han, A. J. Berkley, J. Johansson, P. Bunyk, E. Ladizinsky, S. Govorkov, M. C. Thom, S. Uchaikin, B. Bumble, A. Fung, A. Kaul, A. Kleinsasser, M. H. S. Amin, and D. V. Averin, Probing Noise in Flux Qubits via Macroscopic Resonant Tunneling, Physical Review Letters 101, 117003 (2008).
  * Kairys  _et al._ [2020] P. Kairys, A. D. King, I. Ozfidan, K. Boothby, J. Raymond, A. Banerjee, and T. S. Humble, Simulating the Shastry-Sutherland Ising Model using Quantum Annealing, PRX Quantum 1, 020320 (2020), arXiv:2003.01019 .
  * Crowley  _et al._ [2014] P. J. D. Crowley, T. Đurić, W. Vinci, P. A. Warburton, and A. G. Green, Quantum and classical dynamics in adiabatic computation, Physical Review A 90, 042317 (2014).
  * Albash and Lidar [2015] T. Albash and D. A. Lidar, Decoherence in adiabatic quantum computation, Physical Review A 91, 062320 (2015).
  * Zaborniak and de Sousa [2021] T. Zaborniak and R. de Sousa, Benchmarking Hamiltonian Noise in the D-Wave Quantum Annealer, IEEE Transactions on Quantum Engineering 2, 1 (2021).
  * Lanting  _et al._ [2011] T. Lanting, M. H. S. Amin, M. W. Johnson, F. Altomare, A. J. Berkley, S. Gildert, R. Harris, J. Johansson, P. Bunyk, E. Ladizinsky, E. Tolkacheva, and D. V. Averin, Probing high-frequency noise with macroscopic resonant tunneling, Physical Review B 83, 180502(R) (2011).
  * Dickson  _et al._ [2013] N. G. Dickson, M. W. Johnson, M. H. Amin, R. Harris, F. Altomare, A. J. Berkley, P. Bunyk, J. Cai, E. M. Chapple, P. Chavez, F. Cioata, T. Cirip, P. deBuen, M. Drew-Brook, C. Enderud, S. Gildert, F. Hamze, J. P. Hilton, E. Hoskinson, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Lanting, T. Mahon, R. Neufeld, T. Oh, I. Perminov, C. Petroff, A. Przybysz, C. Rich, P. Spear, A. Tcaciuc, M. C. Thom, E. Tolkacheva, S. Uchaikin, J. Wang, A. B. Wilson, Z. Merali, and G. Rose, Thermally assisted quantum annealing of a 16-qubit problem, Nature Communications 4, 1903 (2013).
  * Marshall  _et al._ [2019] J. Marshall, D. Venturelli, I. Hen, and E. G. Rieffel, Power of Pausing: Advancing Understanding of Thermalization in Experimental Quantum Annealers, Physical Review Applied 11, 044083 (2019).
  * Albash and Marshall [2021] T. Albash and J. Marshall, Comparing relaxation mechanisms in quantum and classical transverse-field annealing, Physical Review Applied 15, 014029 (2021), arXiv:2009.04934 .
  * Chen and Lidar [2020a] H. Chen and D. A. Lidar, Why and when is pausing beneficial in quantum annealing?, arXiv:2005.01888 [quant-ph] (2020a), arXiv:2005.01888 [quant-ph] .
  * Chancellor and Kendon [2021] N. Chancellor and V. Kendon, Experimental test of search range in quantum annealing, Physical Review A 104, 012604 (2021).
  * Shin  _et al._ [2014] S. W. Shin, G. Smith, J. A. Smolin, and U. Vazirani, How "Quantum" is the D-Wave Machine?, arXiv:1401.7087 [quant-ph] (2014), arXiv:1401.7087 [quant-ph] .
  * Muthukrishnan  _et al._ [2016] S. Muthukrishnan, T. Albash, and D. A. Lidar, Tunneling and Speedup in Quantum Optimization for Permutation-Symmetric Problems, Physical Review X 6, 031010 (2016).
  * Albash  _et al._ [2012] T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi, Quantum adiabatic Markovian master equations, New Journal of Physics 14, 123016 (2012).
  * Boixo  _et al._ [2013] S. Boixo, T. Albash, F. M. Spedalieri, N. Chancellor, and D. A. Lidar, Experimental signature of programmable quantum annealing, Nature Communications 4, 2067 (2013).
  * Albash  _et al._ [2015] T. Albash, W. Vinci, A. Mishra, P. A. Warburton, and D. A. Lidar, Consistency Tests of Classical and Quantum Models for a Quantum Annealer, Physical Review A 91, 042314 (2015), arXiv:1403.4228 .
  * Kramers and Wannier [1941] H. A. Kramers and G. H. Wannier, Statistics of the Two-Dimensional Ferromagnet. Part I, Physical Review 60, 252 (1941).
  * Klauder [1979] J. R. Klauder, Path integrals and stationary-phase approximations, Physical Review D 19, 2349 (1979).
  * Kato [1950] T. Kato, On the Adiabatic Theorem of Quantum Mechanics, Journal of the Physical Society of Japan 5, 435 (1950).
  * Chen and Lidar [2020b] H. Chen and D. A. Lidar, HOQST: Hamiltonian Open Quantum System Toolkit, arXiv:2011.14046 [quant-ph] (2020b), arXiv:2011.14046 [quant-ph] .
  * Mishra  _et al._ [2018] A. Mishra, T. Albash, and D. A. Lidar, Finite temperature quantum annealing solving exponentially small gap problem with non-monotonic success probability, Nature Communications 9, 2917 (2018).
  * Zhou  _et al._ [2020] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Quantum Approximate Optimization Algorithm: Performance, Mechanism, and Implementation on Near-Term Devices, Physical Review X 10, 021067 (2020).

##  Appendix A PFC Free-Energy Proof

The transfer matrix representation of the PFC partition function (Eq. (3))
involves singular, non-commuting matrices $\mathbf{W}$ (Eq. (5)) and
$\mathbf{V}=\mathbf{v}^{T}\mathbf{v}$ (where $\mathbf{v}$ is defined in Eq.
(4)), which does not allow for an obvious reduction to an analytical free
energy that is generally defined as

$F=-\lim_{N\rightarrow\infty}\frac{1}{\beta N}\ln\mathcal{Z}\,,$ (23)

where $\mathcal{Z}$ is the partition function. We begin by redefining the
partition function in Eq. (3) to

$\mathcal{Z}=\textrm{Tr}\left(\mathbf{W}^{M-1}\mathbf{V}\right)\,,$ (24)

where $M=N/2$, such that we now take the limit in $M\rightarrow\infty$ due to
$\mathbf{W}$ containing information about the subsystem rather than a single
qubit. Performing an eigen-decomposition on $\mathbf{W}$ to yields a diagonal
matrix of eigenvalues, $\mathbf{D}$, and a matrix of eigenvectors,
$\mathbf{P}$, in the form $\mathbf{W}=\mathbf{P}\mathbf{D}\mathbf{P}^{-1}$.

Figure 9: Probability of being in the ground state for both SVMC-TF (solid
line) and spherical-SVMC (dotted line) as the system scales in size, $M$, for
a PFC with $d=0.1$. Here, the SVMC and spherical-SVMC-TF probabilities are
found from 20,000 samples, which we repeat 50 times and bootstrap to find the
median and $95\%$ confidence intervals for the data point and error bars
respectively. Figure 10: Plots of state probability for being in either the
ground state (solid lines) or any of the $2^{M}$-degenerate $1^{st}$ excited
states (dashed lines) at the end of an anneal. A PFC with system size $M=3$
was evolved using quantum (top row, but now on a log-scale unlike Fig. 7) and
classical (bottom row) dynamics with SVMC-TF and spherical-SVMC. The black
dashed line indicates the probability with random guessing i.e. $1/64$. The
quantum simulations are plotted against anneal time in $ns$, whilst the SVMC
simulations are plotted against the number of sweeps. The closed and open
system dynamics are evolved according to the von-Neumann and adiabatic master
equation respectively (see section III.2). The SVMC-TF and spherical-SVMC
probabilities are found from 20,000 samples, of which we repeat 50 times and
bootstrap to find a median and $95\%$ confidence intervals for the data point
and error bars respectively.

When $\mathbf{W}$ is raised to any power, the decomposition simply becomes
$\mathbf{W}^{M}=\mathbf{P}\mathbf{D}^{M}\mathbf{P}^{-1}$. Given the cyclic
invariance of the trace, the partition function therefore becomes

$\mathcal{Z}=\textrm{Tr}\left(\mathbf{P}\mathbf{D}^{M-1}\mathbf{P}^{-1}\mathbf{V}\right)=\textrm{Tr}\left(\mathbf{D}^{M-1}\mathbf{P}^{-1}\mathbf{V}\mathbf{P}\right)\,.$
(25)

Taking $\mathbf{A}=\mathbf{P}^{-1}\mathbf{V}\mathbf{P}$, and the largest
absolute eigenvalue of $\mathbf{W}$ to be $\lambda_{1}$ (the spectral radius),
the trace summation will yield

$\begin{split}\textrm{Tr}\left(\mathbf{D}^{M-1}\mathbf{A}\right)=\sum_{i=1}^{4}A_{ii}\lambda_{i}^{M-1}=\sum_{i=1}^{4}\left(\sqrt[M-1]{A_{ii}}\lambda_{i}\right)^{M-1}\\\
=\left(\sqrt[M-1]{A_{11}}\lambda_{1}\right)^{M-1}\sum_{i=1}^{4}\left(\frac{\sqrt[M-1]{A_{ii}}\lambda_{i}}{\sqrt[M-1]{A_{11}}\lambda_{1}}\right)^{M-1}\,.\end{split}$
(26)

Given that $\mathbf{W}$ is positive semi-definite, in the limit of
$M\rightarrow\infty$ the trace is simply left with $\lambda_{1}^{M-1}$ as a
non-vanishing term, such that we find the analytical form of the free energy
to be

$F=-\lim_{M\rightarrow\infty}\frac{1}{\beta
M}\ln\textrm{Tr}\left(\mathbf{D}^{M-1}\mathbf{A}\right)=-\frac{1}{\beta}\ln{\lambda_{1}}\,.$
(27)

This eigenvalue can be found symbolically using Pythons SymPy library, such
that

$\begin{split}\lambda_{1}=&e^{2\beta R}\cosh{\beta Rd}+\cosh{\beta R(2-d)}+\\\
&\sqrt{\left(e^{2\beta R}\cosh{\beta Rd}+\cosh{\beta
R(2-d)}\right)^{2}-4\sinh{4\beta R}}\,.\end{split}$ (28)

The validity of this free-energy relies on the fact that our spectral radius
is $\geq 1$, otherwise $\mathbf{W}$ would converge to zero in the limit of
$M\rightarrow\infty$. However, the PFC has bounds of $R>0$ and $0<d<1$, such
that our eigenvalue $\lambda_{1}\geq 4$, and is therefore finite everywhere
except for $\beta R=\infty$.

##  Appendix B Result of SVMC Variants

In this appendix, we present three additional results of the two other SVMC
variants (SVMC-TF and spherical-SVMC) and demonstrate the intermediate effects
in performance when going from the simplest case (SVMC) to the most complex
(spherical-SVMC-TF) variant. The first is of Fig. 9, where we illustrate how
SVMC-TF and spherical-SVMC are affected by increasing the system size,
similarly to what is seen in Fig. 5. Secondly, we see the results for SVMC-TF
and spherical-SVMC (Fig. 10) missing from our comparative results in Fig. 7.
This also includes a log-scale view of the quantum systems tested for
comparison, which highlights what appears to be discontinuities in the
results. These are caused by under-sampling the rapid oscillatory behaviour in
the ground state probability after passing through the minimum gap too quickly
(given that sufficient probability density is placed on the ground state), and
this is amplified for small $d$ (i.e. small minimum gap).

Finally, we show a comparative experiment for the simpler $M=2$ PFC in Fig. 11
that is analogous to the $M=3$ case shown in Fig.7. Here we tune the
perturbative parameter $d$ and show results and simulations for all variants.
As expected, the $M=2$ is easier to solve and we therefore see higher ground
state probabilities in general compared to the $M=3$ case. Furthermore, it
illustrates better how combining both the transverse-field updates and the
extension to the whole Bloch sphere has a enhanced effect and out-performs all
variants for the hardest cases (given the sweep ranges we have tested).

From the variants tested, each individual addition to the SVMC algorithm
brings an improvement to ground state probability for the hardest problem
cases. The helpful addition of the azimuthal component to SVMC can be
explained by the fact that not only does it provide another degree of freedom,
but on average it can reduce the transverse field energy contribution in Eq.
18. This reduces both the size of the minimum gap and moves the position of
the minimum gap to earlier in the anneal, therefore giving SVMC more time
(sweeps) to try and reach the ground state after the minimum gap.

Finally, we would expect the addition of transverse-field updates to provide a
more accurate description of the a quantum annealer (by emulating freeze-out
[26]), but also hinder computation as dynamics would be slow in regimes where
we would need additional dynamics to reach the ground state after the minimum
gap. However, provided that the spin-vector has followed the false minimum and
is near the global minimum of the semi-classical potential (Fig. 3), the spin-
vector would most likely resemble the first excited state closest in Hamming
distance to the ground state after the minimum gap (see Fig. 6). Therefore,
slowed dynamics are preventing the spin-vector from getting lost in the first-
excited-state manifold to some extent, which makes it statistically more
likely for SVMC-TF to reach the ground state.

Figure 11: Plots of state probability for being in either the ground state
(solid lines) or any of the $1^{st}$ excited states (dashed lines) for the
closed and open system quantum simulations and the SVMC variants for a PFC
with system size $M=2$. The black dashed line indicates the probability with
random guessing i.e. $1/16$. The quantum simulations are measured against
anneal time in $ns$, whilst the SVMC simulations are measured against the
number of sweeps used to increment the anneal. Here, the probabilities of the
SVMC variants are found from 20,000 samples, of which we repeat 50 times and
bootstrap to find a median and $95\%$ confidence intervals for the data point
and error bars respectively.


SSfootnotetext: This paper is supported by NSFC 12071352, NSFC 11721101, the
National Key R and D Program of China 2020YFA0713100, and 2022CFB240.

# Liouville Theorem for Harmonic Maps from Riemannian Manifold with Compact
Boundary11footnotemark: 1

Jun Sun<EMAIL_ADDRESS>Xiaobao Zhu<EMAIL_ADDRESS>School of
Mathematics and Statistics, Wuhan University, Wuhan 430072, China  Hubei Key
Laboratory of Computational Science, Wuhan University, Wuhan, 430072, China
School of Mathematics, Renmin University of China, Beijing 100872, China

###### Abstract

In this note we will provide a gradient estimate for harmonic maps from a
complete noncompact Riemannian manifold with compact boundary (which we call
“Kasue manifold”) into a simply connected complete Riemannian manifold with
non-positive sectional curvature. As a consequence, we can obtain a Liouville
theorem. We will also show the nonexistence of positive solutions to some
linear elliptic equations on Kasue manifolds.

###### keywords:

Liouville theorem, gradient estimate, harmonic maps

###### MSC:

[2020]53C20, 53C40

##  1 Introduction

Gradient estimate is a very important technique in geometric analysis and has
attracted much attentions since Yau’s seminal paper ([10]). Yau proved the
following gradient estimate for positive harmonic functions:

######  Theorem 1.1.

([10]) Let $(M,g)$ be an $n$-dimensional, complete Riemmanian manifold
(without boundary). For $K\geq 0$, we assume that $Ric_{M}\geq-(n-1)K$. For
$x_{0}\in M$, let $u:B_{R}(x_{0})\to(0,\infty)$ be a positive harmonic
function. Then we have

$\displaystyle\frac{|\nabla u|}{u}\leq
C_{n}\left(\frac{1}{R}+\sqrt{K}\right),$

on $B_{\frac{R}{2}}(x_{0})$, where $C_{n}$ is a positive constant depending
only on $n$.

A consequence of the gradient estimate is the well-known Liouville theorem,
which states that any positive harmonic function on a complete Riemannian
manifold with nonnegative Ricci curvature is a constant.

Later on, Cheng generalized Yau’s gradient estimate to the harmonic map case
([1]). More precisely, he proved that:

######  Theorem 1.2.

([1]) Let $M$ be a complete Riemannian manifold with Ricci curvature bounded
from below $Ric_{M}\geq-(n-1)K$. Let $N$ be a simply connected complete
Riemannian manifold with non-positive sectional curvature. Let u be a harmonic
map from $M$ to $N$. Assume that $y_{0}\not\in u(B_{R}(0))$. Let $\rho(y)$ be
the distance between $y$ and $y_{0}$ in $N$. Then, if
$b>2\sup\\{\rho(u(x))|x\in B_{R}(0)\\}$, we have

$\displaystyle\sup_{B_{\frac{R}{2}}(0)}\frac{|\nabla
u|(x)}{b^{2}-\rho^{2}(u(x))}\leq\frac{C(1+KR^{2})}{R^{2}b},$

where $C>0$ depends only on $M$ and $N$.

As a consequence of Cheng’s gradient estimate, we can obtain the Liouville
theorem for harmonic maps:

######  Corollary 1.3.

Let $M$ be a complete Riemannian manifold with nonnegative Ricci curvature.
Let $N$ be a simply connected complete Riemannian manifold with non-positive
sectional curvature. Let $u$ be a harmonic map from $M$ to $N$. If the image
of $u$ in $N$ is a bounded set, then it is constant.

Recently, Kunikawa and Sakurai ([6]) derived Yau’s gradient estimate in the
setting that the Riemmannian manifold $M$ is complete noncompact with compact
boundary for harmonic functions with Dirichlet boundary condition (i.e., it is
constant on the boundary). They proved that

######  Theorem 1.4.

([6]) Let $(M,g)$ be an $n$-dimensional, complete Riemmanian manifold with
compact boundary. For $K\geq 0$, we assume that $Ric_{M}\geq-(n-1)K$ and
$H_{\partial M}\geq-(n-1)\sqrt{K}$. Let $u:B_{R}(\partial M)\to(0,\infty)$ be
a positive harmonic function with Dirichlet boundary condition. We assume that
its derivative $u_{\nu}$ in the direction of the outward unit normal vector
$\nu$ is non-negative over $\partial M$. Then we have

$\displaystyle\frac{|\nabla u|}{u}\leq
C_{n}\left(\frac{1}{R}+\sqrt{K}\right),$

on $B_{\frac{R}{2}}(\partial M)$, where $C_{n}$ is a positive constant
depending only on $n$, and $B_{R}(\partial M):=\\{x\in M|d(x,\partial
M)<R\\}$.

In particular, they can obtain the following Liouville theorem:

######  Corollary 1.5.

Let $M$ be a complete Riemannian manifold with compact boundary. We assume
that $Ric_{M}\geq 0$ and $H_{\partial M}\geq 0$. Let $u:M\to(0,\infty)$ be a
positive harmonic function with Dirichlet boundary condition. We assume that
$u_{\nu}\geq 0$ over $\partial M$. Then $u$ is constant.

Actually, a Riemannian manifold $M$ satisfying the assumptions of Corollary
1.5 is classified by Kasue ([5]):

######  Theorem 1.6.

([5]) Let $M$ be a connected, complete noncompact Riemannian manifold with
compact boundary. If $Ric_{M}\geq 0$ and $H_{\partial M}\geq 0$, then
$\partial M$ is connected, and $M$ is isometric to $[0,\infty)\times\partial
M$.

######  Definition 1.1.

We say that $M$ is a Kasue manifold if it satisfies the assumptions in Theorem
1.6.

In this paper we will first prove a Liouville theorem for harmonic maps from a
Kasue manifold to a simply connected Riemannian manifold with non-positive
sectional curvature.

Theorem A: Let $M$ be a Kasue manifold, $N$ be a simply connected complete
Riemannian manifold with non-positive sectional curvature and let $u$ be a
harmonic map from $M$ to $N$ with Dirichlet boundary condition (i.e., $u$ is
constant on the boundary). If the image of $u$ in $N$ is a bounded set, then
it is constant.

The Liouville theorem will follow from the general gradient estimate (see
Theorem 2.4). The idea of the proof of the gradient estimate follows from that
of Cheng ([1]) which applied the maximum principle to appropriately chosen
test function. The main difference is to deal with the case that the maximum
of the test function is achieved on the boundary of $M$. For this reason, our
gradient estimate holds with specially chosen $y_{0}$, while Cheng’s gradient
estimate holds for any chosen $y_{0}$. However, this is enough to guarantee
the validity of the Liouville theorem.

The other aim of this paper is to show the nonexistence of positive solution
to some elliptic equations on Kasue manifolds. Precisely,

Theorem B: Let $(M,g)$ be a Kasue manifold and suppose $h\in C^{2}(M)$
satisfies

$\displaystyle\Delta h\geq 0,~{}~{}0\leq h\not\equiv
0~{}~{}\text{and}~{}~{}h_{\nu}|_{\partial M}\leq 0.$

Then the equation

$\displaystyle\begin{cases}&\Delta
u+hu=0~{}~{}~{}~{}\text{in}~{}~{}~{}~{}M,\\\ &u_{\nu}|_{\partial M}\geq
0\end{cases}$ (1.1)

does not admit a positive solution.

To end the introduction, we want to say some words on the motivation of our
results. When we read the article of Kunikawa and Sakurai [6], we noticed that
the Laplacian comparison theorem also holds for the distance function from the
boundary. Relating to our knowledge on gradient estimates, we realized that we
can derived some Liouville properties for elliptic equations and harmonic
maps. Then our results appeared. The method of our proof is standard. We
believe that one can use it to study other elliptic equations and parabolic
equations on Kasue manifolds. Here, we refer the reader to [3] for gradient
estimates for a nonlinear parabolic equation on Kasue manifolds.

##  2 Proof of Theorem A

In this section, we will prove Theorem A. We will use the classical idea to
prove the gradient estimate (Theorem 2.4). To that purpose, we first recall
the Laplacian comparison theorem for the distance function from a compact
hypersurface in the Riemannian manifold $M$.

Let $(M,g)$ be an $n$-dimensional complete noncompact Riemannian manifold with
compact boundary $\partial M$. The distance function from the boundary
$r_{\partial M}:M\to{\mathbb{R}}$ is defined as

$\displaystyle r_{\partial M}:=d(\cdot,\partial M),$

which is smooth outside of the cut locus for the boundary Cut$\partial M$
([9]).

For $K,\Lambda\in{\mathbb{R}}$, we denote by $s_{K,\Lambda}(t)$ the unique
solution to the Jacobi equation $\varphi^{\prime\prime}(t)+K\varphi(t)=0$ with
initial conditions $\varphi(0)=1$ and $\varphi^{\prime}(0)=-\Lambda$. We have
the following Laplacian comparison theorem:

######  Theorem 2.1.

([4]) For $K,\Lambda\in{\mathbb{R}}$, we assume $Ric_{M}\geq(n-1)K$ and
$H_{\partial M}\geq(n-1)\Lambda$. Then we have

$\displaystyle\Delta r_{\partial
M}\leq(n-1)\frac{s^{\prime}_{K,\Lambda}(r_{\partial
M})}{s_{K,\Lambda}(r_{\partial M})}$

outside of Cut$\partial M$.

The following Bochner formula will be used in the derivation of gradient
estimate:

######  Lemma 2.2.

([2]) Let $(M^{n},g)$ and $(N^{m},h)$ be two Riemannian manifolds and $u:M\to
N$ be a harmonic map, then the energy density of $u$ satisfies the following
formula:

$\frac{1}{2}\Delta|\nabla u|^{2}=|\nabla du|^{2}+\langle Ric^{M}\nabla
u,\nabla u\rangle-\langle
Rm^{N}(u_{\alpha},u_{\beta})u_{\alpha},u_{\beta}\rangle,$

where $Ric^{M}$ and $Rm^{N}$ are the Ricci curvature of $(M,g)$ and Riemannian
curvature of $(N,h)$, respectively.

We will also need the following Reilly type formula which was first proved for
functions by Reilly ([8]). For the map case, we refer to [7].

######  Lemma 2.3.

([7]) Let $(M^{n},g)$ be an $n$-dimensional Riemannian manifold with compact
boundary $\partial M$, and $(N^{m},h)$ be a Riemannian manifold. Let $\nu$ be
the unit outer normal vector of $\partial M$ in $M$. Then for all smooth map
$u:(M,g)\to(N,h)$, we have

$\displaystyle\left(|\nabla u|^{2}\right)_{\nu}=$ $\displaystyle 2\langle
du(\nu),\tau^{\partial
M}(u)-\tau(u)\rangle-2H|du(\nu)|^{2}+2\langle\nabla(du(\nu)),du\rangle_{\partial
M}$ $\displaystyle-2\langle du\circ A,du\rangle_{\partial M},$

where $Av:=\nabla^{M}_{v}\nu$ for $v\in T\partial M$ is the Weingarten
operator, $H$ is the mean curvature of $\partial M$ in $M$ with respect to
$\nu$, and $\tau^{\partial M}(u)$ and $\tau(u)$ are the tension fields of
$u_{\partial M}:\partial M\to N$ and $u:M\to N$, respectively.

Now we can prove the gradient estimate.

######  Theorem 2.4.

Let $(M,g)$ be an $n$-dimensional, complete Riemmanian manifold with compact
boundary and $N$ be a simply connected complete Riemannian manifold with non-
positive sectional curvature. For $K\geq 0$, we assume that
$Ric_{M}\geq-(n-1)K$ and $H_{\partial M}\geq-(n-1)\sqrt{K}$. Let $u:M\to N$ be
a harmonic map with Dirichlet boundary condition. Assume that the image of $u$
is a bounded set in $N$. Then we can choose $y_{0}\not\in u(M)$ so that, if we
let $\rho(y)$ be the distance between $y$ and $y_{0}$ in $N$, then for any
$R>0$ we have

$\displaystyle\sup_{B_{\frac{R}{2}}(\partial M)}\frac{|\nabla
u|(x)}{b^{2}-\rho^{2}(u(x))}\leq\frac{C(1+\sqrt{K}R)}{bR},$

for some constant $b>2\sup\\{\rho(u(x))|x\in M\\}$, where $C>0$ depends only
on $M$ and $N$.

###### Proof.

Let $(M^{n},g)$ be an $n$-dimensional complete noncompact Riemannian manifold
with compact boundary $\partial M$, and $(N^{m},h)$ be a simply connected
Riemannian manifold with non-positive sectional curvature. Let $u:M\to N$ be a
smooth harmonic map with Dirichlet boundary condition. We assume that the
image of $u$, denoted by $u(M)$, is bounded in $N$, and $y_{1}:=u(\partial
M)\in u(M)$. In the following, the constant $C$ will denote a constant
depending only on $n$, which may vary from line to line.

Now we fix a point $y_{0}\not\in u(M)$, which will be specified later. Let
$\rho(y)$ be the distance between $y$ and $y_{0}$ in $N$. Fix
$b>2\sup\\{\rho(u(x))|x\in M\\}$. As in [1], we define

$\displaystyle\phi(x)=\frac{|\nabla u(x)|^{2}}{(b^{2}-\rho^{2}(u(x)))^{2}}.$

Then we have

$\displaystyle\nabla\phi(x)=\frac{\nabla(|\nabla
u(x)|^{2})}{(b^{2}-\rho^{2}(u(x)))^{2}}+2\frac{|\nabla
u(x)|^{2}\nabla(\rho^{2}(u(x)))}{(b^{2}-\rho^{2}(u(x)))^{3}},$ (2.1)

and

$\displaystyle\Delta\phi(x)=$ $\displaystyle\frac{\Delta(|\nabla
u(x)|^{2})}{(b^{2}-\rho^{2}(u(x)))^{2}}+4\frac{\nabla(|\nabla
u(x)|^{2})\cdot\nabla(\rho^{2}(u(x)))}{(b^{2}-\rho^{2}(u(x)))^{3}}$
$\displaystyle+2\frac{|\nabla
u(x)|^{2}\Delta(\rho^{2}(u(x)))}{(b^{2}-\rho^{2}(u(x)))^{3}}+6\frac{|\nabla
u(x)|^{2}|\nabla(\rho^{2}(u(x)))|^{2}}{(b^{2}-\rho^{2}(u(x)))^{4}}.$ (2.2)

The Bochner formula for harmonic maps together with our assumptions on
curvatures implies that

$\displaystyle\Delta|\nabla u|^{2}\geq 2|\nabla du|^{2}-2(n-1)K|\nabla
u|^{2}.$

By standard argument using Cauchy-Schwartz inequality and Hessian comparison
theorem, we can obtain the following inequality:

$\displaystyle\Delta\phi(x)\geq\frac{4|\nabla
u(x)|^{4}}{(b^{2}-\rho^{2}(u(x)))^{3}}-\frac{2(n-1)K|\nabla
u(x)|^{2}}{(b^{2}-\rho^{2}(u(x)))^{2}}+\frac{2\nabla\phi\cdot\nabla\rho^{2}}{b^{2}-\rho^{2}(u(x))}.$
(2.3)

Now for any $R>0$, we define a function $F:B_{R}({\partial M})\to{\mathbb{R}}$
by

$\displaystyle F(x):=(R^{2}-r^{2}_{\partial M}(x))^{2}\phi(x).$

We assume that $F$ achieves its maximum at some point $x_{0}\in B_{R}(\partial
M)$.

Case 1: $x_{0}\in B_{R}(\partial M)\backslash\partial M$. In this case, using
Calabi’s trick, we may assume that $x_{0}$ does not belong to Cut$\partial M$.
Therefore, at $x_{0}$, it holds that

$\displaystyle\nabla F(x_{0})=0$

and

$\displaystyle\Delta F(x_{0})\leq 0.$

Hence we have at $x_{0}$ that

$\displaystyle\frac{\nabla\phi}{\phi}=\frac{4r_{\partial M}\nabla r_{\partial
M}}{R^{2}-r_{\partial M}^{2}}$

and

$\displaystyle\frac{\Delta\phi}{\phi}-\frac{8r_{\partial M}\nabla r_{\partial
M}\cdot\nabla\phi}{(R^{2}-r_{\partial M}^{2})\phi}-\frac{2\Delta
r^{2}_{\partial M}}{R^{2}-r_{\partial M}^{2}}+\frac{8r^{2}_{\partial
M}}{(R^{2}-r_{\partial M}^{2})^{2}}\leq 0.$

It follows that

$\displaystyle\frac{\Delta\phi}{\phi}-\frac{24r^{2}_{\partial
M}}{(R^{2}-r_{\partial M}^{2})^{2}}-\frac{2\Delta r^{2}_{\partial
M}}{R^{2}-r_{\partial M}^{2}}\leq 0.$

Notice that $s_{-K,-\sqrt{K}}(t)=e^{\sqrt{K}t}$. The Laplacian comparison
theorem (Theorem 2.1) and our assumptions imply that

$\displaystyle\Delta r^{2}_{\partial M}=2r_{\partial M}\Delta r_{\partial
M}+2\leq 2(n-1)\sqrt{K}r_{\partial M}+2\leq C(1+\sqrt{K}r_{\partial M}).$

Combining with (2.3), we obtain that

$\displaystyle 0\geq$
$\displaystyle\frac{\Delta\phi}{\phi}-\frac{24r^{2}_{\partial
M}}{(R^{2}-r_{\partial M}^{2})^{2}}-\frac{2\Delta r^{2}_{\partial
M}}{R^{2}-r_{\partial M}^{2}}$ $\displaystyle\geq$ $\displaystyle
4(b^{2}-\rho^{2})\phi-2(n-1)K+\frac{8r_{\partial M}\nabla r_{\partial
M}\cdot\nabla\rho^{2}}{(R^{2}-r_{\partial M}^{2})(b^{2}-\rho^{2})}$
$\displaystyle-\frac{24r^{2}_{\partial M}}{(R^{2}-r_{\partial
M}^{2})^{2}}-\frac{C(1+\sqrt{K}r_{\partial M})}{R^{2}-r_{\partial M}^{2}}.$

Because

$\displaystyle\rho(u(x))\leq\frac{b}{2},$

and

$\displaystyle|\nabla(\rho^{2}(u(x)))|\leq 2\rho|\nabla^{N}\rho||\nabla u|\leq
b|\nabla u|,$

we have

$\displaystyle 3b^{2}\phi-\frac{8r_{\partial M}b|\nabla u|}{(R^{2}-r_{\partial
M}^{2})(b^{2}-\rho^{2})}-\frac{24r^{2}_{\partial M}}{(R^{2}-r_{\partial
M}^{2})^{2}}-\frac{C(1+\sqrt{K}r_{\partial M})}{R^{2}-r_{\partial
M}^{2}}-2(n-1)K\leq 0.$

Multiplying the above inequality through $(R^{2}-r_{\partial M}^{2})^{2}$, we
have

$\displaystyle 3b^{2}F-8r_{\partial M}bF^{\frac{1}{2}}-C(1+\sqrt{K}r_{\partial
M})R^{2}-2(n-1)KR^{4}\leq 0,$

which yields that

$\displaystyle\sup_{B_{\frac{R}{2}}(\partial M)}F^{\frac{1}{2}}(x)\leq
F^{\frac{1}{2}}(x_{0})\leq\frac{CR(1+\sqrt{K}R)}{b}.$

In particular, we have

$\displaystyle\sup_{B_{\frac{R}{2}}(\partial M)}\frac{|\nabla
u|}{R^{2}-\rho^{2}}\leq F^{\frac{1}{2}}(x_{0})\leq\frac{C(1+\sqrt{K}R)}{bR}.$

Case 2: $x_{0}\in\partial M$. In this case, we know that

$\displaystyle F_{\nu}(x_{0})\geq 0.$

Since $r_{\partial M}(x_{0})=0$, we can easily see that

$\displaystyle\phi_{\nu}(x_{0})=\frac{F_{\nu}(x_{0})}{R^{4}}\geq 0.$ (2.4)

Choose local orthonormal frame $\\{e_{1},\cdots,e_{n-1},\nu\\}$ of $M$ along
$\partial M$ around $x_{0}$ so that $T\partial M$ is spanned by
$\\{e_{1},\cdots,e_{n-1}\\}$. By our assumption, $u(\partial M)=y_{1}$. This
implies that at $x_{0}$, we have

$\displaystyle|du(\nu)|=|\nabla u|,\ \ du(e_{i})=0,\ for\ 1\leq i\leq n-1.$

Also we have $\tau(u)=0$ since $u$ is harmonic. Hence, Lemma 2.2 implies that

$\displaystyle(|\nabla u|^{2})_{\nu}(x_{0})=-2H(x_{0})|\nabla u|^{2}(x_{0}).$

Therefore, (2.4) gives us at $x_{0}$ that

$\displaystyle 0\leq$ $\displaystyle\phi_{\nu}=\frac{(|\nabla
u|^{2})_{\nu}}{(b^{2}-\rho^{2})^{2}}+\frac{2|\nabla u|^{2}(\rho^{2}\circ
u)_{\nu}}{(b^{2}-\rho^{2})^{3}}$ $\displaystyle=$
$\displaystyle-2H\frac{|\nabla
u|^{2}}{(b^{2}-\rho^{2})^{2}}+\frac{4\rho|\nabla
u|^{2}\langle\nabla^{N}\rho,du(\nu)\rangle}{(b^{2}-\rho^{2})^{3}}.$ (2.5)

Notice that $\nabla^{N}\rho(u(x_{0}))=\nabla^{N}\rho(y_{1})$ is the radial
direction from $y_{0}$ to $y_{1}$. Since $du(\nu)$ is independent of the
choice of $y_{0}$, we can choose $y_{0}\not\in u(M)$ so that $\nabla\rho$ is
in the direction of $-du(\nu)$. Since $u(M)$ is bounded in $N$, $y_{0}$ and
$b$ can also be chosen so that

$\displaystyle\rho(y_{1})\geq\inf\\{\rho(u(x))|x\in
M\\}\geq\frac{2}{3}\sup\\{\rho(u(x))|x\in M\\}\geq\frac{1}{4}b.$

With this choice of $y_{0}$, we see that
$\langle\nabla^{N}\rho,du(\nu)\rangle=-|du(\nu)|=-|\nabla u|$. Therefore, we
have from (2) that

$\displaystyle b\phi^{\frac{3}{2}}(x_{0})\leq-2H(x_{0})\phi(x_{0})\leq
2(n-1)\sqrt{K}\phi(x_{0})$

so that

$\displaystyle
F^{\frac{1}{2}}(x_{0})=R^{2}\phi^{\frac{1}{2}}(x_{0})\leq\frac{C\sqrt{K}R^{2}}{b}\leq\frac{CR(1+\sqrt{K}R)}{b}.$

In particular, we have

$\displaystyle\sup_{B_{\frac{R}{2}}(\partial M)}\frac{|\nabla
u|}{R^{2}-\rho^{2}}\leq F^{\frac{1}{2}}(x_{0})\leq\frac{C(1+\sqrt{K}R)}{bR}.$

This finishes the proof of the theorem. ∎

###### Proof of Theorem A.

The main theorem follows by setting $K=0$ and letting $R$ tend to infinity in
Theorem 2.4. ∎

##  3 Proof of Theorem B

In this section, we shall give the proof of Theorem B, which is also based on
a gradient estimate.

###### Proof of Theorem B.

We first assume that $(M,g)$ is an $n$-dimensional, complete Riemmanian
manifold with compact boundary, $Ric_{M}\geq-(n-1)K$ and $H_{\partial
M}\geq-(n-1)\sqrt{K}$ for some $K\geq 0$.

Suppose $u$ is a positive solution to equation

$\displaystyle\Delta u+hu=0~{}~{}\text{on}~{}~{}M.$

Let

$\displaystyle w=\frac{|\nabla u|^{2}}{u^{2}}.$

By choosing a local orthonormal system, we calculate the equation of $w$ as
follows:

$\displaystyle w_{j}=\frac{2u_{i}u_{ij}}{u^{2}}-\frac{2u_{i}^{2}u_{j}}{u^{3}}$
(3.1)

and

$\displaystyle\Delta
w=w_{jj}=\frac{2u_{ij}^{2}}{u^{2}}+\frac{2u_{i}u_{ijj}}{u^{2}}-\frac{4u_{i}u_{ij}u_{j}}{u^{3}}-\frac{4u_{i}u_{ij}u_{j}}{u^{3}}-\frac{2u_{i}^{2}u_{jj}}{u^{3}}+\frac{6u_{i}^{2}u_{j}^{2}}{u^{4}}.$

Then by the Ricci formula and (3.1), we have

$\displaystyle\Delta w=$ $\displaystyle
2\left(\frac{u_{ij}}{u}-\frac{u_{i}u_{j}}{u^{2}}\right)^{2}-2\nabla
w\cdot\nabla\log u+\frac{2u_{i}(\Delta
u)_{i}}{u^{2}}+\frac{2R_{ij}u_{i}u_{j}}{u^{2}}-\frac{2u_{i}^{2}\Delta
u}{u^{3}}$ $\displaystyle\geq$ $\displaystyle\frac{2}{n}\left(\frac{\Delta
u}{u}-\frac{|\nabla u|^{2}}{u^{2}}\right)^{2}-2\nabla w\cdot\nabla\log
u-2\nabla h\cdot\nabla\log u-2(n-1)Kw$ $\displaystyle=$
$\displaystyle\frac{2}{n}(w+h)^{2}-2\nabla(w+h)\cdot\nabla\log u-2(n-1)Kw.$

Since we assume $\Delta h\geq 0$, one obtains

$\displaystyle\Delta(w+h)\geq\frac{2}{n}(w+h)^{2}-2\nabla(w+h)\cdot\nabla\log
u-2(n-1)Kw.$ (3.2)

Let $\psi$ be a smooth cut-off function supported in $B_{R}(\partial M)$,
satisfying the following properties:

(1) $\psi(x)=\psi(r_{\partial M}(x))$ and $\psi^{\prime}\leq 0$;

(2) $0\leq\psi\leq 1$ and $\psi\equiv 1$ in $B_{R/2}(\partial M)$;

(3) $\frac{|\psi^{\prime}|}{\psi^{a}}\leq\frac{C_{a}}{R}$,
$\frac{|\psi^{\prime\prime}|}{\psi^{a}}\leq\frac{C_{a}}{R^{2}}$ when
$a\in(0,1)$.

Case 1: If $(w+h)\psi$ attains its maximum at some point $x_{1}\in
B_{R}(\partial M)\setminus\partial M$, by using Calabi’s argument we can
assume w.l.o.g. that $x_{1}\not\in\text{Cut}\partial M$, then

$\nabla[(w+h)\psi](x_{1})=0,~{}~{}~{}~{}~{}~{}\Delta[(w+h)\psi](x_{1})\leq 0.$

Calculating directly and using (3.2), one has at $x_{1}$ that

$\displaystyle 0\geq$
$\displaystyle\psi\Delta(w+h)+2\nabla(w+h)\cdot\nabla\psi+(w+h)\Delta\psi$
$\displaystyle\geq$
$\displaystyle\frac{2}{n}(w+h)^{2}\psi+2(w+h)\nabla\psi\cdot\nabla\log
u-2(n-1)Kw\psi$
$\displaystyle-2(w+h)\frac{|\nabla\psi|^{2}}{\psi}+(w+h)\Delta\psi.$ (3.3)

We multiply inequality (3) by $\psi(x_{1})$ and estimate the new right-hand-
side term by term.

Firstly, since $h$ is nonnegative, it follows from Cauchy’s inequality and the
properties of $\psi$ that

$\displaystyle 2(w+h)\psi\nabla\psi\cdot\nabla\log u\geq$
$\displaystyle-(w+h)^{3/2}\psi|\psi^{\prime}|$ $\displaystyle\geq$
$\displaystyle-\frac{1}{4n}(w+h)^{2}\psi^{2}-C(n)\left(\frac{|\psi^{\prime}|}{\psi^{1/2}}\right)^{4}$
$\displaystyle\geq$
$\displaystyle-\frac{1}{4n}(w+h)^{2}\psi^{2}-C(n)\frac{1}{R^{4}}.$ (3.4)

Secondly, we have

$\displaystyle-2(n-1)Kw\psi^{2}\geq-\frac{1}{4n}(w+h)^{2}\psi^{2}-C(n)K^{2}.$
(3.5)

Thirdly, one has

$\displaystyle-2(w+h)|\nabla\psi|^{2}\geq$
$\displaystyle-\frac{1}{4n}(w+h)^{2}\psi^{2}-C(n)\left(\frac{|\psi^{\prime}|^{2}}{\psi}\right)^{2}$
$\displaystyle\geq$
$\displaystyle-\frac{1}{4n}(w+h)^{2}\psi^{2}-C(n)\frac{1}{R^{4}}.$ (3.6)

Lastly, by the Laplace comparison theorem (Theorem 2.1) we have

$\displaystyle(w+h)\psi\Delta\psi\geq$
$\displaystyle-\frac{1}{4n}(w+h)^{2}\psi^{2}-C(n)(\Delta\psi)^{2}$
$\displaystyle=$
$\displaystyle-\frac{1}{4n}(w+h)^{2}\psi^{2}-C(n)[\psi^{\prime\prime}+\psi^{\prime}\Delta
r_{\partial M}]^{2}$ $\displaystyle\geq$
$\displaystyle-\frac{1}{4n}(w+h)^{2}\psi^{2}-C(n)[|\psi^{\prime\prime}|+|\psi^{\prime}|\frac{n-1}{r_{\partial
M}}(1+\sqrt{K})r_{\partial M}]^{2}$ $\displaystyle\geq$
$\displaystyle-\frac{1}{4n}(w+h)^{2}\psi^{2}-C(n)\left(\frac{1}{R^{4}}+K^{2}\right).$
(3.7)

Substituting (3)-(3) into (3), one obtains

$\displaystyle\frac{1}{n}[(w+h)\psi]^{2}(x_{1})\leq
C(n)\left(\frac{1}{R^{4}}+K^{2}\right).$

Therefore,

$\displaystyle[(w+h)\psi](x)\leq[(w+h)\psi](x_{1})\leq
C(n)\left(\frac{1}{R^{2}}+K\right).$

Since $\psi\equiv 1$ in $B_{R/2}(\partial M)$, we have

$\displaystyle\sup_{B_{R/2}(\partial M)}w\leq\sup_{B_{R/2}(\partial
M)}(w+h)\leq C(n)\left(\frac{1}{R^{2}}+K\right).$

Case 2: If $(w+h)\psi$ attains its maximum at some point $x_{2}\in\partial M$,
one has at $x_{2}$ that

$[(w+h)\psi]_{\nu}\geq 0,$

and hence

$w_{\nu}\geq(w+h)_{\nu}=\frac{1}{\psi}\left([(w+h)\psi]_{\nu}-(w+h)\psi_{\nu}\right)=[(w+h)\psi]_{\nu}\geq
0.$

The Dirichlet boundary condition for $u$ and the assumption $u_{\nu}\geq 0$
tell us that $|\nabla u|=u_{\nu}$ on $\partial M$. It follows by the Reilly
formula (see [8] or Proposition 2.3 in [6]) that

$\displaystyle 0\leq$ $\displaystyle w_{\nu}\leq(|\nabla\log
u|^{2})_{\nu}=2(\log u)_{\nu}(\Delta u-(\log u)_{\nu}H)$ $\displaystyle=$
$\displaystyle 2\frac{u_{\nu}}{u}\left(\frac{\Delta u}{u}-\frac{|\nabla
u|^{2}}{u^{2}}-\frac{u_{\nu}}{u}H\right)$ $\displaystyle=$ $\displaystyle
2\frac{u_{\nu}}{u}\left(-h-\frac{|\nabla
u|^{2}}{u^{2}}-\frac{u_{\nu}}{u}H\right)$ $\displaystyle\leq$ $\displaystyle
2\frac{u_{\nu}}{u}\left(-w-\frac{u_{\nu}}{u}H\right).$

Therefore,

$\displaystyle w\leq\max\\{0,-w^{1/2}H\\}\leq(n-1)\sqrt{K}w^{1/2}.$

We then have

$w\leq C(n)K.$

Combining the above two cases, we obtain finally that

$\displaystyle\sup_{B_{R/2}(\partial M)}\frac{|\nabla
u|^{2}}{u^{2}}=\sup_{B_{R/2}(\partial M)}w\leq
C(n)\left(\frac{1}{R^{2}}+K\right).$ (3.8)

When $M$ is a Kasue manifold, we have $K=0$. Then by passing $R$ tend to
infinity in (3.8), one has $u$ is a constant. Recalling that $h$ is
nonnegative and not identical to $0$, the equation (1.1) can not admit a
positive constant solution. So we arrive at that (1.1) has no positive
solution and complete the proof of Theorem B. ∎

## References

  * Cheng [1980] S.-Y. Cheng, Liouville theorem for harmonic maps, Geometry of the Laplace operator, Proc. Sympos. Pure Math. 36, Amer. Math. Soc., Providence, R.I., 1980, 147–151.
  * Eells and Sampson [1964] Jr. J. Eells and J. H. Sampson, Harmonic mappings of Riemannian manifolds, Amer. J. Math. 86 (1964), 109–160.
  * Fu and Wu [2022] X. Fu and J.-Y. Wu, Gradient estimates for a nonlinear parabolic equation with Dirichlet boundary condition, Kodai Math. J. 45 (2022), 96–109.
  * Kasue [1982] A. Kasue, A Laplacian comparison theorem and function theoretic properties of a complete Riemannian manifold, Japan. J. Math. (N.S.) 8 (1982), 309–341.
  * Kasue [1983] A. Kasue, Ricci curvature, geodesics and some geometric properties of Riemannian manifolds with boundary, J. Math. Soc. Japan 35 (1983), 117–131.
  * Kunikawa and Sakurai [2008] K. Kunikawa and Y. Sakurai, Yau and Souplet-Zhang type gradient estimates on Riemannian manifolds with boundary under Dirichlet boundary condition, Proc. Amer. Math. Soc. 150 (2022), 1767–1777.
  * Matei [2000] A. Matei, Gap phemomena for $p$-harmonic maps, Ann. Global Anal. Geom. 18 (2000), 541–554.
  * Reilly [1977] R. Reilly, Applications of the Hessian operator in a Riemannian manifold, Indiana Univ. Math. J. 26 (1977), 459–472.
  * Sakurai [2017] Y. Sakurai, Rigidity of manifolds with boundary under a lower Ricci curvature bound, Osaka J. Math. 54 (2017), 85–119.
  * Yau [1975] S.-T. Yau, Harmonic functions on complete Riemannian manifolds, Comm. Pure Appl. Math. 28 (1975), 201–228.