text:
# Dynamics of quantum adiabatic evolution algorithm for Number Partitioning

## I  Introduction

Since the discovery by Shor Shor  nearly a decade ago of a quantum algorithm for efficient integer factorization there has been a rapidly growing interest in the development of new quantum algorithms capable of solving computational problems that are practically intractable on classical computers. Perhaps the most notable example is that of a combinatorial optimization problem (COP). In the simplest case the task in COP is to minimize the cost function (“energy”) $`E_𝐳`$ defined on a set of $`2^n`$ binary strings $`𝐳=\{z_1,\mathrm{},z_n\}`$ $`z_j=0,1`$, each containing $`n`$ bits. In quantum computation this cost function corresponds to a Hamiltonian $`H_P`$

$`H_P={\displaystyle \underset{𝐳}{}}E_𝐳|𝐳𝐳|`$ (1)
$`|𝐳=|z_1_1|z_2_2\mathrm{}|z_n_n.`$
where $`z_j=0,1`$ and the summation is over $`2^n`$ states $`|𝐳`$ forming the computational basis of a quantum computer with $`n`$ qubits. State $`|z_j_j`$ of the $`j`$-th qubit is an eigenstate of the Pauli matrix $`\widehat{\sigma }_z`$ with eigenvalue $`S_j=12z_j(S_j=\pm 1)`$. It is clear from the above that the ground state of $`H_P`$ encodes the solution to the COP with cost function $`E_𝐳`$.

COPs have a direct analogy in physics, related to finding ground states of classical spin glass models. In the example above bits $`z_j`$ correspond to Ising spins $`S_j`$. The connection between the properties of frustrated disordered systems and the structure of the solution space of complex COPs has been noted first by Fu and Anderson Anderson . It has been recognized Parizi  that many of the spin glass models are in almost one-to-one correspondence with a number of COPs from theoretical computer science that form the so-called NP-complete class Garey . This class contains hundreds of the most common computationally hard problems encountered in practice, such as constraint satisfaction, traveling salesmen, and integer programming. NP-complete problems are characterized in the worst cases by exponential scaling of the running time or memory requirements with the problem size $`n`$. A special property of the class is that any NP-complete problem can be converted into any other NP-complete problem in polynomial time on a classical computer; therefore, it is sufficient to find a deterministic algorithm that can be guaranteed to solve all instances of just one of the NP-complete problems within a polynomial time bound. It is widely believed, however, that such an algorithm does not exist on a classical computer; whether it exists on a quantum computer is one of the central open questions. Ultimately, one can expect that the behavior of new quantum algorithms for COPs and their complexity will be closely related to the properties of quantum spin glasses.

Recently, Farhi and co-workers suggested a new quantum algorithm for solving combinatorial optimization problems which is based on the properties of quantum adiabatic evolution Farhi . Running of the algorithm for several NP-complete problems has been simulated on a classical computer using a large number of randomly generated problem instances that are believed to be computationally hard for classical algorithms FarhiSc ; FarhiSat ; FarhiCli ; Hogg . Results of these numerical simulations for relatively small size of the problem instances ( $`n`$ 20) suggest a quadratic scaling law of the run time of the quantum adiabatic algorithm with $`n`$. Furthermore, it was shown in Vazirani02  that the previous query complexity argument that led to the exponential lower bound for unstructured search Bennett  cannot be used to rule out the polynomial time solution of NP-complete Satisfiability problem by the quantum adiabatic algorithm.

In Vazirani01 ; annealing ; Vazirani02 ; farhimeas ; BS  special symmetric cases of COP were considered where symmetry of the problem allowed the authors to describe the true asymptotic behavior ($`n\mathrm{}`$) of the algorithm. In certain examples considered in Farhi ; annealing  the quantum adiabatic algorithm finds the solution in time polynomial in $`n`$ while simulated annealing requires exponential time. This effect occurs due to the special connectivity properties of the optimization problems that lead to the relatively large matrix elements for the spin tunneling in transverse magnetic field between different valleys during the quantum adiabatic algorithm. In the examples considered in annealing  the tunneling matrix element scales polynomially with $`n`$. On the other hand, in simulated annealing different valleys are connected via classical activation processes for spins with probabilities that scale exponentially with $`n`$. It was also shown for certain simplified examples farhimeas ; BS , that quantum adiabatic algorithm can be modified to completely suppress the tunneling barriers even if the corresponding classical cost function has local minima well separated in the space of spin configurations.

However, so far there are no study on the true asymptotic behavior of the algorithm for the general case of randomly generated hard instances of NP-complete problems. Also there are no analysis of the limitations of the quantum adiabatic computation arising from the intrinsic properties of disorder and frustration in this problems. Such analysis is of the central interest in this paper.

In Sec. II we introduce the random Number Partitioning problem and describes conditional cost distributions (neighborhood properties) in this problem. In Sec. III we describe the concept of quantum adiabatic computation applied to combinatorial optimization problems and introduce a Green function method for the analysis of the minimum gap. In Sec. IV we describe the effect of quantum diffusion in the algorithm dynamics, derive the scaling for the minimum gap and the complexity of the algorithm for the random Number Partitioning problem. We also obtain the scaling of the minimum gap numerically from the form of the cumulative density of the adiabatic eigenvalues at the avoided-crossing point. In Sec. V we discuss the results of the simulations of the time-dependent Schrödinger equation to simulate quantum adiabatic computation for Number Partitioning and obtain its complexity numerically.

## II  Number Partitioning Problem

Number Partitioning Problem (NPP) is one of the six basic NP-complete problems that are at the heart of the theory of NP-completeness Garey . It can be formulated as a combinatorial optimization problem: Given a sequence of positive numbers $`\{a_1,\mathrm{},a_n\}`$ find a partition, i.e. two disjoint subsets $`𝒜`$ and $`𝒜^{}`$, such that the residue

$$E=\left|\underset{a_j𝒜}{}a_j\underset{a_j𝒜^{}}{}a_j\right|$$
(2)
is minimized. In NPP we search for the bit strings $`𝐳=\{z_1,\mathrm{},z_n\}`$ (or corresponding Ising spin configurations $`𝐒=\{S_1,\mathrm{},S_n\}`$) that minimize the energy or cost function $`E_𝐳`$

$$E_𝐳=\left|\mathrm{\Omega }_𝐒\right|,\mathrm{\Omega }_𝐒=\underset{j=1}{\overset{n}{}}a_jS_j,S_j=12z_j,$$
(3)
where $`S_j=1`$ ($`z_j=0`$) if $`a_j𝒜`$ and $`S_j=1`$ ($`z_j=1`$) if $`a_j𝒜^{}`$. The partition $`𝐒`$ with minimum residue can also be viewed as the ground state of the Ising spin glass, $`\mathrm{\Omega }_𝐒^2`$, corresponding to the Mattis-like antiferromagnetic coupling, $`J_{ij}=a_ia_j`$.

NPP has many practical applications including multiprocessor scheduling ms , cryptography cr , and others. The best deterministic heuristical algorithm for NPP, the differencing method of Karmakar and Karp Karp , can find with high probability solutions whose energies are of the order $`1/n^{\alpha \mathrm{log}n}`$ for some $`\alpha >0`$. The interest in NPP also stems from the remarkable failure of a standard simulated annealing algorithm for the energy function (3) to find good solutions, as compared with the solutions found by deterministic heuristics johnson . The apparent reason for this failure is due to the existence of order $`2^n`$ local minima whose energies are of the order of $`1/n`$ Ferreira  which undermines the usual strategy of exploring the space of the spin configurations $`𝐒`$ through single spin flips.

The computational complexity of random instances of NPP depends on the number of bits $`b`$ needed to encode the numbers $`a_j`$. In what follows we will analyze NPP with independent, identically distributed (i.i.d.) random b-bit numbers $`a_j`$. Numerical simulations show Gent ; Korf ; MertensCompleteAnytime  that solution time grows exponentially with $`n`$ for $`nb`$ then decreases steeply for $`nb`$ (phenomenon of “peaking”) and eventually grows polynomially for $`nb`$. The transition from the “hard” to computationally “easy” phases at $`nb`$ has features somewhat similar to phase transitions in physical systems Mertens98 . The detailed theory of the phase transition in NPP was given in Refs. Borgs1 ; Borgs2 . If one keeps the parameter $`\xi =b/n`$ fixed and lets $`n\mathrm{}`$ then instances of NPP corresponding to $`\xi >`$ 1 will have no perfect partitions with high probability. On the other hand for $`\xi <`$ 1 number of perfect partitions will grow exponentially with $`n`$. Transitions of this kind were observed in various NP-complete problems AI96 . In what follows we will focus on the computationally hard regime $`\xi 1`$.

### II.1  Distribution of signed partition residues

The values of individual energies are random and depend on the particular instance of NPP (i.e., the set of numbers $`a_j`$). However on a coarse-grained scale (i.e. after averaging over individual energy separations) the form of the typical energy distribution is described by some universal function for randomly generated problem instances. We introduce for a given set of randomly sampled numbers $`a_j`$ a coarse-grained distribution function of signed partition residues $`\mathrm{\Omega }_𝐳`$ (3)

$$P(\mathrm{\Omega })=2^n\frac{1}{\mathrm{\Delta }\mathrm{\Omega }}_{\mathrm{\Omega }\mathrm{\Delta }\mathrm{\Omega }/2}^{\mathrm{\Omega }+\mathrm{\Delta }\mathrm{\Omega }/2}𝑑\eta \underset{𝐳\{0,1\}^n}{}\delta (\eta \mathrm{\Omega }_𝐳).$$
(4)
Here $`\delta (x)`$ is the Dirac delta-function; the sum is over $`2^n`$ bit-strings $`𝐳`$ and $`2^n`$ is a normalization factor. In (4) we average over an interval $`\mathrm{\Delta }\mathrm{\Omega }`$ of the partition residues whose size is chosen self-consistently, $`\mathrm{\Delta }\mathrm{\Omega }2^n/P(\mathrm{\Omega })`$. Using (3) we can rewrite (4) in the form

$`P(\mathrm{\Omega })={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{\mathrm{}}}𝑑s\zeta \left({\displaystyle \frac{\mathrm{\Delta }\mathrm{\Omega }s}{2}}\right)I(s)\mathrm{cos}(\mathrm{\Omega }s),`$ (5)
$`I(s)={\displaystyle \underset{j=1}{\overset{n}{}}}\mathrm{cos}(a_js),\zeta (x)=\mathrm{sin}(x)/x.`$
Here $`\zeta (x)`$ is a window function that imposes a cut-off in the integral (5) at $`s2/\mathrm{\Delta }\mathrm{\Omega }`$. For large $`n`$ this integral can be evaluated using the steepest descent method. In the following we shall assume that the $`b`$-bit numbers $`a_j`$ are distributed inside of the unit interval $`[0,1]`$ and are integer multiples of $`2^b`$, the smallest number that can be represented with available number of bits $`b`$. We note that for large $`n`$ the function $`I(s)`$ has sharp maxima (minima) with width $`n^{1/2}`$ at the points $`s_k=k\pi \mathrm{\hspace{0.17em}2}^b,k=0,1,\mathrm{};`$ $`|I(s_k)|=1`$. Only one saddle point at $`s=0`$ contributes to the integral in (5) due to coarse-graining of the distribution (4). Indeed, it will be seen below that the window size $`2/\mathrm{\Delta }\mathrm{\Omega }`$ can be chosen to obey the conditions $`1n^{1/2}/\mathrm{\Delta }\mathrm{\Omega }2^n`$. Therefore in the case of high-precision numbers, $`bn`$, saddle-points $`s_k`$ with $`k>0`$ lie far outside the window and their contributions can be neglected (see also Appendix A). On the other hand the window function $`\zeta (x)`$ can be replaced by unity while computing the contribution from the saddle-point at $`s=0`$. Finally we obtain for $`|\mathrm{\Omega }|n`$ (cf. Mertens2000 )

$`P(\mathrm{\Omega })={\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2(0)n}}}\mathrm{exp}\left({\displaystyle \frac{\mathrm{\Omega }^2}{2\sigma ^2(0)n}}\right)+𝒪(n^{3/2})`$
$`\sigma ^2(0)={\displaystyle \frac{1}{n}}{\displaystyle \underset{j=1}{\overset{n}{}}}a_j^2(En).`$ (6)
The coarse-grained distribution $`P(\mathrm{\Omega })`$ depends on the set of $`a_j`$’s through a single self-averaging quantity $`\sigma (0)`$ (cf. Mertens98 ).

One can also introduce the distribution $`\stackrel{~}{P}(E)`$ of cost values (energies) $`E_𝐳=|\mathrm{\Omega }_𝐳|`$. Due to the obvious symmetry of the NPP, the cost function $`E_𝐳`$ in (3) does not change after flipping signs of all spins, $`S_jS_j`$. Therefore

$$\stackrel{~}{P}(E)=1/2P(\pm E).$$
(7)
We emphasize that, according to Eq. (6) for a typical set of high-precision numbers $`a_j`$ the energy spectrum in NPP is quasi-continuous, and there are only two scales present in the distribution $`\stackrel{~}{P}(E)`$: one is a “microscopic” scale given by the characteristic separation of the individual partition energies, $`E_{\mathrm{min}}`$, and another is given by the mean partition energy $`E`$ (or the distribution width $`E^2^{1/2}`$)

$$E_{\mathrm{min}}\sigma (0)n^{1/2}\mathrm{\hspace{0.17em}\hspace{0.17em}2}^n,E^2=\frac{\pi }{2}E^2=n\sigma ^2(0).$$
(8)
This justifies the choice for $`\mathrm{\Delta }\mathrm{\Omega }`$ above that corresponds to coarse-graining over many individual energy level separations.

We note that the distribution $`P(\mathrm{\Omega })`$ (6) is Gaussian for $`En`$ and can be understood in terms of a random walk with coordinate $`\mathrm{\Omega }`$ using Eq. (3). The walk begins at the origin, $`\mathrm{\Omega }=0`$, and makes a total of $`n`$ steps. At the $`j`$-th step $`\mathrm{\Omega }`$ moves to the right or to the left by “distance” $`2a_j`$ if $`S_j=1`$ or $`S_j=1`$, respectively. In the asymptotic limit of large $`n`$ the result (6) corresponds to equal probabilities of right and left moves and the distribution of step lengths coinciding with that of the set of numbers $`\{2a_j\}`$.

Finally, the energy distribution function $`P(E)`$ of the form (6),(7) was previously obtained by Mertens Mertens2000  using explicit averaging over the random instances of NPP. He also computed the partition function $`Z(T)`$ for a given instance of NPP at a small finite temperature $`T`$ using the steepest-descent method and summation over the saddle-points $`s_k=k\pi \mathrm{\hspace{0.17em}2}^b`$ similar to our discussion above Mertens98  (in his analysis $`k_BT`$ played a role similar to our regularization factor $`\mathrm{\Delta }\mathrm{\Omega }`$ in (4),(5)).

We emphasize however, that the approach in Ref. Mertens98  based on $`Z(T)`$ is necessarily restricted to the analysis of the “static” properties of NPP at $`E2^n`$, i.e., the phase transition in the number of perfect partitions Mertens98  when the control parameter $`\xi =n/b`$ crosses a critical value. On the other hand distribution $`P(\mathrm{\Omega })`$ (4) introduces at finite energies, as well as the conditional distribution introduced in the next section also allow us to directly study the intrinsic dynamical properties of the problem in question such as the dynamics of its quantum optimization algorithms.

### II.2 Conditional distribution of signed partition residues

Consider the set of bit-strings $`𝐳^{}`$ obtained from a given string $`𝐳`$ by flipping $`r`$ bits. The conditional distribution of the partition residues $`\mathrm{\Omega }_𝐳^{}`$ (3) in the $`r`$-neighborhood of $`𝐳`$ can be characterized by its moments:

$$\mathrm{\Omega }^k=\left(\genfrac{}{}{0pt}{}{n}{r}\right)^1\underset{𝐳^{}\{0,1\}^n}{}\left(\mathrm{\Omega }_𝐳^{}\right)^k\delta _{r,D(𝐳^{},𝐳)},k=1,2,\mathrm{}$$
(9)
Here $`\delta _{m,l}`$ is a Kronecker delta and function $`D(𝐳,𝐳^{})`$ computes the number of bits that take different values in the bit-strings $`𝐳`$ and $`𝐳^{}`$. It is the so-called Hamming distance between the two strings

$$D(𝐳,𝐳^{})=\underset{j=1}{\overset{n}{}}\left|z_jz_j^{}\right|.$$
(10)
The Hamming distance $`r=D(𝐳,𝐳^{})`$ between the bit-strings is directly related to the overlap factor $`q`$ between the corresponding spin configurations often used in the theory of spin glasses Parizi ; Mertens2000 :

$$q=\frac{1}{n}\underset{j=1}{\overset{n}{}}S_jS_j^{}=1\frac{2}{n}D(𝐳,𝐳^{}).$$
(11)
(in what follows we shall use both quantities $`r`$ and $`q`$). For $`k=1,2`$ in (9) one obtains after straightforward calculation the first and second moments of the conditional distribution

$`\mathrm{\Omega }=q\mathrm{\Omega }_𝐳,`$ (12)
$`\mathrm{\Omega }^2\mathrm{\Omega }^2=n\sigma ^2(q)\left(1+{\displaystyle \frac{1}{n1}}\right)\left(1{\displaystyle \frac{1}{n}}{\displaystyle \frac{\mathrm{\Omega }_𝐳^2}{E^2}}\right),`$ (13)
$`\sigma (q)=\sigma (0)(1q^2)^{1/2},q1{\displaystyle \frac{2r}{n}},`$ (14)
where $`\sigma (0)`$ and $`E^2`$ are given in (6) and (8), respectively.

The conditional distribution of $`\mathrm{\Omega }_𝐳^{}`$ can also be defined in a way similar to (4)

$$P_{r,𝐳}(\mathrm{\Omega }^{})=\left(\genfrac{}{}{0pt}{}{n}{r}\right)^1\frac{1}{\mathrm{\Delta }\mathrm{\Omega }^{}}_{\mathrm{\Omega }^{}\mathrm{\Delta }\mathrm{\Omega }^{}/2}^{\mathrm{\Omega }^{}+\mathrm{\Delta }\mathrm{\Omega }^{}/2}𝑑\eta \underset{𝐳^{}\{\mathrm{𝟎},\mathrm{𝟏}\}^𝐧}{}\delta (\eta \mathrm{\Omega }_𝐳^{})\delta _{r,D(𝐳^{},𝐳)}$$
(15)
where averaging is over the small interval $`\mathrm{\Delta }\mathrm{\Omega }^{}`$ that, however, includes many individual values of $`\mathrm{\Omega }_z^{}`$ for a given $`r`$. It is clear from (12),(13) that the first two moments of $`P_{r,𝐳}(\mathrm{\Omega }^{})`$ depend on $`𝐳`$ only via the value of $`\mathrm{\Omega }_𝐳`$. This does not hold true, however, for the higher-order moments that depend on other functions of $`𝐳`$ as well. For example, $`\mathrm{\Omega }^3`$ involves the quantity $`_{j=1}^na_j^3(12z_j)`$, etc.

Our main observation is that in the asymptotic limit of large $`n`$ the conditional distribution $`P_{r,𝐳}(\mathrm{\Omega }^{})`$ is well-described by the first two moments (12),(13). Then, according to the discussion above, its dependence on $`𝐳`$ is only via $`\mathrm{\Omega }_𝐳`$. The detailed study of the higher moments (9) will be done elsewhere. Here we use the following intuitive approach relevant for analysis of the computational complexity of the quantum adiabatic algorithm for the NPP. We average $`P_{r,𝐳}(\mathrm{\Omega }^{})`$ over the strings $`𝐳`$ with residues $`\mathrm{\Omega }_𝐳`$ inside a small interval $`\mathrm{\Delta }\mathrm{\Omega }`$ (containing, however, many levels $`\mathrm{\Omega }_𝐳`$). After such averaging the result, $`P_r(\mathrm{\Omega }^{}|\mathrm{\Omega })`$, can be written in the form

$`P_r(\mathrm{\Omega }^{}|\mathrm{\Omega })={\displaystyle \frac{P_r(\mathrm{\Omega }^{},\mathrm{\Omega })}{P(\mathrm{\Omega })}},`$ (16)
$`P_r(\mathrm{\Omega }^{},\mathrm{\Omega })=2^n{\displaystyle \frac{1}{\mathrm{\Delta }\mathrm{\Omega }}}{\displaystyle _{\mathrm{\Omega }\mathrm{\Delta }\mathrm{\Omega }/2}^{\mathrm{\Omega }+\mathrm{\Delta }\mathrm{\Omega }/2}}𝑑\eta {\displaystyle \underset{𝐳\{0,1\}^n}{}}\delta (\eta \mathrm{\Omega }_𝐳)P_{r,𝐳}(\mathrm{\Omega }^{}),`$ (17)
where $`P(\mathrm{\Omega })`$ is given in (6). We note that Eq. (16) formally coincides with the Bayesian rule expressing the conditional distribution function $`P_r(\mathrm{\Omega }^{}|\mathrm{\Omega })`$ through the 2-point (joint) distribution function $`P_r(\mathrm{\Omega }^{},\mathrm{\Omega })`$ and the single-point distribution $`P(\mathrm{\Omega })`$ (6).

The explicit form of $`P_r(\mathrm{\Omega }|\mathrm{\Omega })`$ is derived in Appendix B in a manner similar to the derivation of $`P(E)`$ in Sec. II.1. The results are presented in Eqs. (72) and (73). They show that $`P_r(\mathrm{\Omega }^{},\mathrm{\Omega })`$ in the limit $`n1`$ is indeed well described by its first two moments that corresponds precisely to the expressions given in Eqs. (12),(13) above. From this we conclude that

$$P_{r,𝐳}(\mathrm{\Omega }^{})=P_r(\mathrm{\Omega }^{}|\mathrm{\Omega }_𝐳).$$
(18)
In the case $`r=1`$ there are $`n`$ strings $`𝐳^{}`$ at a Hamming distance 1 from the string $`𝐳`$. Partition energies corresponding to these strings equal $`|\mathrm{\Omega }_𝐒2a_jS_j|,\mathrm{\hspace{0.17em}1}jn`$ (cf. (3)). After the coarse-graining over the energy scale $`𝒪(1/n)`$ in the range, $`|\mathrm{\Omega }|,|\mathrm{\Omega }^{}|n`$, the conditional distribution $`P_{r,𝐳}`$ is a step function in the interval $`\mathrm{\Omega }_𝐳\mathrm{\Omega }^{}[2,2]`$. For $`r=n1`$ one has the same form of the distribution but for $`\mathrm{\Omega }_𝐳+\mathrm{\Omega }^{}`$. Both results correspond to nearly equal distribution of spins between between $`\pm 1`$ values. Then in the range of energies $`|\mathrm{\Omega }^{}|,|\mathrm{\Omega }_𝐳|1`$ one has:

$$P_{r,𝐳}(\mathrm{\Omega }^{})\overline{P}_r=1/2+𝒪\left(\frac{1}{n}\right),r=1,n1(n1)$$
(19)
For $`r,nr1`$ distribution $`P_{r,𝐳}(\mathrm{\Omega }^{})`$ has a Gaussian form with a broad maximum at $`\mathrm{\Omega }^{}=q\mathrm{\Omega }_𝐳`$ (cf. Eqs. (12),(13),(72)) . Near the maximum we have:

$$P_{r,𝐳}(\mathrm{\Omega }^{})\overline{P}_r=\frac{1}{\sqrt{2\pi n\sigma ^2(q)}},|\mathrm{\Omega }^{}|,|\mathrm{\Omega }_𝐳|n^{1/2}\sigma (q).$$
(20)
We studied the conditional distribution in NPP numerically as well (see Fig.1 and Sec.B). The results are in good agreement with theory even for modest values of $`n30`$.

The characteristic spacing between the values of the partition residues in the subset of strings $`𝐳^{}`$ with $`D(𝐳^{},𝐳)=r`$ is $`1/(\overline{P}_r\left(\genfrac{}{}{0pt}{}{n}{r}\right))`$ for not too large $`E_𝐳,E_𝐳^{}`$ (see above). This spacing decreases exponentially with the magnitude of the string overlap factor, $`|q|=|(n2r)/n|`$. The hierarchy of the subsets corresponding to different values of $`|q|`$ form a specific structure of NPP. We note that the distribution of partition residues within the hierarchy is nearly independent of the ancestor string $`𝐳`$ in a broad range of energies $`E^{}n^{1/2}`$ where $`P_{r,𝐳}(E^{})\overline{P}_r`$. One can see that the magnitude of the overlap factor $`q`$ between two strings with energies within a given interval $`[0,E]`$ is limited by some typical value $`\overline{q}`$ satisfying the following equation:

$$E\left(\genfrac{}{}{0pt}{}{n}{r}\right)\overline{P}_r=1,|\overline{q}|=12\frac{r}{n}.$$
(21)
The smaller $`E`$ is, the smaller $`|\overline{q}|`$ is: strings that are close in energy are far away in the configuration space. This property gives rise to an exponentially large number of local minima for small values of $`E_𝐳`$ that are far apart in the configuration space. For example, strings with $`E_𝐳E_{\mathrm{min}}`$ typically correspond to $`|q|=𝒪(1/n)`$, they can be obtained from each other only by simultaneously flipping clusters with $`n/2`$ spins.

Eq. (21) describes the dynamics of a local search heuristic (e.g., simulated annealing). It shows that the average cost value $`E`$ during the search decreases no faster than $`𝒪(1/M)`$ where $`M=𝒪\left(\left(\genfrac{}{}{0pt}{}{n}{r}\right)\right)`$ is the number of generated configurations. This result coincides with that obtained in Mertens2000  using a different approach. It says that any classical local search heuristic for NPP cannot be faster than random search. Indeed, during local search the information about the “current” string $`𝐳`$ with $`E_𝐳1`$ is being lost, on average, after one spin flip (cf. Eqs. (19),(20)). We show below that precisely this property of NPP also leads to the complexity of the quantum adiabatic algorithm corresponding to that of a quantum random search.

We note that one can trivially break the symmetry of NPP mentioned above by introducing an extra number $`a_0`$ and placing it, say, in the subset $`𝒜`$. In this case different partition energies will still be encoded by spin configurations $`𝐒=\{S_1,\mathrm{},S_n\}`$ (or corresponding bit-strings $`𝐳`$) with $`\mathrm{\Omega }_𝐒=a_0+_{j=1}^nS_ja_j`$ and $`E_𝐳=|\mathrm{\Omega }_𝐒|`$ (cf. 3). We shall adopt this approach in the analysis of the performance of the quantum adiabatic algorithm for NPP given below.

## III  Quantum Adiabatic Evolution Algorithm

In the quantum adiabatic algorithm Farhi  one specifies the time-dependent Hamiltonian $`H(t)=\stackrel{~}{H}(t/T)`$

$$\stackrel{~}{H}(\tau )=(1\tau )V+\tau H_P,$$
(22)
where $`\tau =t/T`$ is dimensionless “time”. This Hamiltonian guides the quantum evolution of the state vector $`|\psi (t)`$ according to the Schrödinger equation $`i|\psi (t)t=H(t)|\psi (t)`$ from $`t=0`$ to $`t=T`$, the run time of the algorithm (we let $`\mathrm{}=1`$). $`H_P`$ is the “problem” Hamiltonian given in (1). $`V`$ is a “driver” Hamiltonian, that is designed to cause transitions between the eigenstates of $`H_P`$. In this algorithm one prepares the initial state of the system $`\psi (0)`$ to be the ground state of $`\stackrel{~}{H}(0)=V`$. In the simplest case

$$V=\underset{j=1}{\overset{n}{}}\sigma _x^j,|\psi (0)=2^{n/2}\underset{𝐳}{}|𝐳,$$
(23)
where $`\sigma _x^j`$ is a Pauli matrix for $`j`$-th qubit. Consider instantaneous eigenstates $`|\varphi _\eta (\tau )`$ of $`\stackrel{~}{H}(\tau )`$ with energies $`\lambda _\eta (\tau )`$ arranged in nondecreasing order at any value of $`\tau (0,1)`$

$$\stackrel{~}{H}|\varphi _\eta =\lambda _\eta |\varphi _\eta ,\eta =0,1,\mathrm{},2^n1.$$
(24)
Provided the value of $`T`$ is large enough and there is a finite gap for all $`t(0,T)`$ between the ground and excited state energies, $`g(\tau )=\lambda _1(\tau )\lambda _0(\tau )>0`$, quantum evolution is adiabatic and the state of the system $`|\psi (t)`$ stays close to an instantaneous ground state, $`|\varphi _0(t/T)`$ (up to a phase factor). Because $`H(T)=H_P`$ the final state $`|\psi (T)`$ is close to the ground state $`|\varphi _0(\tau =1)`$ of the problem Hamiltonian. Therefore a measurement performed on the quantum computer at $`t=T(\tau =1)`$ will find one of the solutions of COP with large probability.

There is a broad class of COPs from theoretical Computer Science where the number of distinct values of a cost function scales polynomially in the size of an input $`n`$. An example is the Satisfiability problem in which the cost $`E_𝐳`$ of a given string $`𝐳`$ equals the number of constrains violated by the string. For those problems, the spectrum of $`H(\tau )`$, at the beginning ($`\tau 0`$) and at the end ($`\tau 1`$) of the algorithm, consists of a polynomial number of well-separated energy levels. Quantum transitions away from the adiabatic ground state occur most likely near the avoided-crossing points $`\tau \tau ^{}`$ where the energy gap $`g(\tau )`$ reaches its minima Hogg . Near the avoided-crossing points, the spectrum of $`H(\tau )`$ is quasi-continuous, with the separation between individuals eigenvalues scaled down with $`n`$. The probability of a quantum transition, $`1|\psi (t)|\varphi _0(t/T)|_{t=T}^2`$, is small provided that

$$T\frac{|\varphi _1|\stackrel{~}{H}_\tau |\varphi _0|_{\tau =\tau ^{}}}{g_{\mathrm{min}}^2},g_{\mathrm{min}}=\underset{0\tau 1}{\mathrm{min}}\left[\lambda _1(\tau )\lambda _0(\tau )\right],$$
(25)
($`\stackrel{~}{H}_\tau d\stackrel{~}{H}/d\tau `$). The fraction in (25) gives an estimate for the required runtime of the algorithm and the task is to find its asymptotic behavior in the limit of large $`n1`$. The numerator in (25) is less than the largest eigenvalue of $`\stackrel{~}{H}_\tau =H_PV`$, typically polynomial in $`n`$ Farhi . However, $`g_{\mathrm{min}}`$ can scale down exponentially with $`n`$ and in such cases the runtime of the quantum adiabatic algorithm will grow exponentially with the size of COP.

### III.1  Implementation of QAA for NPP

As suggested in Farhi  the quantum adiabatic algorithm can be recast within the conventional quantum computing paradigm using the technique introduced by Lloyd Lloyd . Continuous-time quantum evolution can be approximated by a time-ordered product of unitary operators, $`e^{i(1\tau _k)V\delta }e^{i\tau _kH_P\delta }`$, corresponding to small time intervals $`(t_k,t_k+\delta )`$. Operator $`e^{i(1\tau _k)V\delta }`$ typically corresponds to a sequence of 1- or 2-qubit gates (cf. (23)). Operator $`e^{i\tau _kH_P\delta }`$ is diagonal in the computational basis $`|𝐳`$ and corresponds to phase rotations by angles $`E_𝐳\delta `$. Since in the case $`nb`$, the average separation between the neighboring values of $`E_𝐳`$ is $`1/P(E)=𝒪(2^n)`$, the quantum device would need to support a very high precision in its physical parameters (like external fields, etc.) to control small $`𝒪(2^n)`$ differences in phases. Since this precision scales with $`n`$ exponentially it would strongly restrict the size of an instance of NPP that could be solved on such a quantum computer. This technical restriction is generic for COPs that involve a quasi-continuous spectrum of cost-function values. Among the other examples are many Ising spin glass models in physics (e.g., the Sherrington-Kirkpatrick model Parizi ). To avoid this restriction we introduce a new oracle-type cost function $`_𝐳`$ that returns a set of values

$$_𝐳=c(\mathrm{\Omega }_𝐳),c(x)\{\epsilon _0,\epsilon _1,\mathrm{},\epsilon _M\}(\epsilon _{k+1}>\epsilon _k),$$
(26)
that can be stored using a relatively small number of bits $`𝒪(\mathrm{log}n)`$. For example, we can divide an interval of partition energies $`(0,B)`$, $`B=_{j=0}^{j=n}a_j`$ into bins whose sizes grow exponentially with the energy. Then the new cost will take one value per bin

$`c(x)=\epsilon _kM+k\mathrm{for}\omega _k|x|<\omega _{k+1},`$
$`\omega _k=(2^k1)\mathrm{\Delta },k=0,\mathrm{},M.`$ (27)
The last bin is $`\omega _M|\mathrm{\Omega }_𝐳|B`$ where we have $`_𝐳=\epsilon _M=0`$. The value of the cutoff $`\omega _MB`$ is discussed below. In this example the Hilbert space of $`2^n`$ states $`|𝐳`$ is divided into $`M+1`$ subspaces $`_k`$, each determined by Eq. (27) for a given $`k`$

$$H_P=\underset{k=0}{\overset{M}{}}\epsilon _k\underset{𝐳_k}{}|𝐳𝐳|.$$
(28)
Note that subspace $`_0`$ contains the solution(s) to NPP. Dimension $`d_0`$ of $`_0`$ is controlled by the value of $`\mathrm{\Delta }`$ in (27) which is another control parameter of the algorithm. We set $`\mathrm{\Delta }=2^nK/P(0)`$ where the integer $`Kd_01`$ is independent of $`n`$ and determines how many times on average one needs to repeat the quantum algorithm in order to obtain the solution to NPP with probability close to 1.

Operator $`H_P`$ projects any state $`|\psi `$ onto the states with partition residues in the range $`0|\mathrm{\Omega }_𝐳|<\omega _M`$. If we choose

$$1\omega _ME,$$
(29)
then the distribution function (6) is nearly uniform for $`|\mathrm{\Omega }_𝐳|\omega _M`$. Therefore the dimensions of the subspaces $`_k`$ grow exponentially with $`k`$: $`d_k=d_0\mathrm{\hspace{0.17em}2}^k`$ for $`k<M`$. This simplification would not affect the complexity of a quantum algorithm that spends most of its time in “annealing” the system to much smaller partition residues, $`\omega _M|\mathrm{\Omega }_𝐳|E_{\mathrm{min}}=𝒪(n^{1/2}\mathrm{\hspace{0.17em}2}^n)`$.

We note that the new discrete-valued cost function defined in (27) is non-local. Unlike problems such as Satisfiability, it cannot be represented by a sum of terms each involving a small number of bits. To implement a unitary operator $`e^{i\tau _kH_P\delta }`$ with $`H_P`$ given in (28) one needs to implement the following classical function on a quantum computer

$$_𝐳=\mathrm{\Theta }(\omega _M|\mathrm{\Omega }_𝐳|)\left[\mathrm{log}_2\left(\frac{\mathrm{\Delta }+|\mathrm{\Omega }_𝐳|}{\mathrm{\Delta }+\omega _M}\right)\right],\mathrm{\Omega }_𝐳=\underset{j=1}{\overset{n}{}}a_j(12z_j).$$
(30)
Here $`[x]`$ denotes the integer part of a number $`x`$; $`\mathrm{\Theta }(x)`$ is the theta-function ($`\mathrm{\Theta }(x)=1`$ for $`x0`$ and $`\mathrm{\Theta }(x)=0`$ for $`x<0`$). The implementation of (30) with quantum circuits involves, among other things, the addition of $`n`$ numbers together with their signs to compute $`\mathrm{\Omega }_𝐳`$, and taking the discrete logarithm of a $`b`$-bit number with respect to base $`2`$. These operations can be performed using a number of quantum gates that is only polynomial in $`n`$ and $`b`$ (cf. Shor  for the implementation of the discrete logarithm).

Since the implementation of a cost function (26),(30) does not add an exponential overhead to the complexity of QAA the feasibility of this algorithm for NPP depends on the scaling of the minimum gap $`g_{\mathrm{min}}`$ with $`n`$.

### III.2  Stationary Schrödinger equation for adiabatic eigenstates

We now solve the stationary Schrödinger equation (24) and obtain the minimum gap $`g_{\mathrm{min}}`$ (25) in the asymptotic limit $`n\mathrm{}`$. To proceed we need to introduce a new basis of states $`|𝐱=|x_1_1|x_2_2\mathrm{}|x_n_n`$ where state $`|x_j_j`$ is an eigenstate of the Pauli matrix $`\widehat{\sigma }_x`$ for the $`j`$-th qubit with eigenvalue $`12x_j=\pm 1`$. Driver Hamiltonian $`V`$ can be written in the following form:

$$V=\underset{m=0}{\overset{n}{}}V_m^m,^m=\underset{x_1+\mathrm{}+x_n=m}{}|𝐱𝐱|.$$
(31)
For a particular case given in Eq. (23) we have $`V_m=2mn`$. Matrix elements of $`^m`$ in a basis of states $`|𝐳`$ depend only on the Hamming distance $`D(𝐳,𝐳^{})`$ between the strings $`𝐳`$ and $`𝐳^{}`$

$$𝐳|^m|𝐳^{}=I_{D(𝐳,𝐳^{})}^m,$$
(32)
$$I_r^m=2^n\underset{q=0}{\overset{nr}{}}\underset{p=0}{\overset{r}{}}\left(\genfrac{}{}{0pt}{}{nr}{q}\right)\left(\genfrac{}{}{0pt}{}{r}{p}\right)(1)^p\mathrm{\Delta }_{m,q+p}.$$
(33)
We now rewrite Eq. (24) in the form

$$|\varphi =\frac{\tau }{\lambda \alpha V}H_P|\varphi ,\alpha \alpha (\tau )=1\tau ,$$
(34)
(we drop the subscript $`\eta `$ indicating the number of a quantum state and also the argument $`\tau `$ in $`\varphi `$ and $`\lambda `$). From (27)-(34) we obtain the equation for the amplitudes $`\varphi _𝐳=𝐳|\varphi `$ in terms of the coefficients $`I_r^m`$

$$\left[1\tau G_0c(\mathrm{\Omega }_𝐳)\right]\varphi _𝐳=\frac{\tau \mathrm{\Phi }\mathrm{\hspace{0.17em}2}^n}{\lambda \alpha V_0}+\tau \underset{𝐳^{}𝐳}{}G_{D(𝐳,𝐳^{})}\varphi _𝐳^{}c(\mathrm{\Omega }_𝐳),$$
(35)
$$\mathrm{\Phi }=\underset{𝐳^{}}{}c(\mathrm{\Omega }_𝐳^{})\varphi _𝐳^{},$$

$$G_rG_r(\lambda )=\underset{m=1}{\overset{n}{}}\frac{I_r^m}{\lambda \alpha V_m},0rn.$$

Here we separated out a “symmetric” term $`2^n\mathrm{\Phi }`$ corresponding to the coupling between the states $`|𝐳`$ via the projection operator $`^0`$ (31).

## IV  Minimum gap analysis

### IV.1  Coarse-graining of the transition matrix

We now make a key observation that $`\varphi _𝐳`$ in (35) can be determined based on the properties of the conditional distribution $`P_{r,𝐳}(E)`$ (15) and the form of the Green function $`G_r(\lambda )`$. We sum the Green function $`G_{D(𝐳,𝐳^{})}`$ over all possible transitions from a given state $`𝐳^{}`$ to states $`𝐳^{}𝐳`$ with energy $`\epsilon _k`$. For not too large partition residues of the initial and final states we obtain

$$\underset{𝐳_k,𝐳𝐳^{}}{}G_{D(𝐳,𝐳^{})}(\lambda )F_k(\lambda )+f_{𝐳^{},k}(\lambda )$$
(36)
$$F_k(\lambda )=\frac{\mu s(\lambda )}{2^{Mk}},s(\lambda )=_0^n𝑑r\frac{\sigma (0)}{\sigma \left(12r/n\right)}\left(\genfrac{}{}{0pt}{}{n}{r}\right)G_r(\lambda )$$
(37)
$$|\mathrm{\Omega }_𝐳^{}|,|\mathrm{\Omega }_𝐳|E,\mu =\frac{2\omega _M}{\pi E}.$$
(38)
Function $`\sigma (q)`$ above is defined in (14) and $`f_{𝐳^{},k}(\lambda )`$ is a small correction described below. In function $`s(\lambda )`$ we replaced summation over the integer values of $`r`$ by an integral. It can be evaluated using the explicit form of $`G_r(\lambda )`$ that decays rapidly with $`r`$. In what follows we will be interested in the region $`|\lambda \alpha V_0|1`$ where

$$2\alpha G_r(\lambda )=\left(\genfrac{}{}{0pt}{}{n}{r}\right)^1\underset{m=1}{\overset{nr}{}}\frac{2^n\left(\genfrac{}{}{0pt}{}{n}{m+r}\right)}{m}2^n\left(\mathrm{ln}r+\gamma \right).$$
(39)
($`\gamma `$ is Euler’s constant) and $`s(\lambda )\mathrm{ln}2/(2\alpha )`$. We note that

$$2\alpha G_r(\lambda )\left[(n/2r)\left(\genfrac{}{}{0pt}{}{n}{r}\right)\right]^1,n/2r1.$$
(40)
Therefore the integrand in $`s(\lambda )`$ is a smooth function of $`r`$ for $`rn/2`$ and quickly decays to zero for $`rn/2`$. The contribution to the integral in $`s(\lambda )`$ from the range of $`rn`$ is small ($`𝒪((r/n)^{1/2}`$).

We note that term $`F_k`$ in (36) provides an “entropic” contribution to the sum in (36). It comes from the large number of states $`𝐳_k`$ corresponding to large Hamming distances $`r`$ from the state $`𝐳^{}`$, $`1rn/2`$. Each state contributes a small weight, $`G_r\left(\genfrac{}{}{0pt}{}{n}{r}\right)^1`$, and number of states for a given $`r`$ is large, $`(\omega _{k+1}\omega _k)\left(\genfrac{}{}{0pt}{}{n}{r}\right)\overline{P}_r1`$. Here $`(\omega _{k+1}\omega _k)`$ is an energy bin for the subspace $`_k`$ and $`\overline{P}_r`$ is the conditional density of states described in Sec. II. The size of the bin scales down exponentially with $`k`$ (cf. (27)) and so does the entropic term $`F_k`$. Below a certain cross-over value of $`k`$ one has $`|F_k||f_{𝐳^{},k}(\lambda )|`$. In this case the dominant contribution to the sum (36) comes from the states $`𝐳`$ with small $`r=D(𝐳,𝐳^{})1`$. In particular for $`k=0`$ one can obtain

$$f_{𝐳^{},0}(\lambda )G_1(\lambda )\underset{𝐰_0}{}\delta _{1,D(𝐳^{},𝐰)}+𝒪(n^3),$$
(41)
where the higher-order terms correspond to $`D(𝐳^{},𝐰)`$ 2. According to (39), $`|G_1(\lambda )|n^2`$ and therefore $`|f_{𝐳,0}|`$ is exponentially larger than the entropic term, $`|F_0|\omega _0d_0\mathrm{\hspace{0.17em}2}^n`$. We note that, unlike the entropic term, $`f_{𝐳^{},0}`$ strongly depends on $`𝐳^{}`$ due to the discreteness of the partition energy spectrum ($`\omega _0n1`$). E.g., depending on a state $`𝐳^{}`$, in this case there could be either one or none of the states $`𝐰_0`$ in the sum (41) satisfying $`D(𝐳^{},𝐰)=1`$.

### IV.2  Extended and localized eigenstates

Based on the discussion above we look for solution of Eq. (35) in the following form:

$$\varphi _𝐳=v(\mathrm{\Omega }_𝐳)+u_𝐳,𝐳_0,$$
(42)
where we have explicitly separated out a part of the wavefunction $`v(\mathrm{\Omega }_𝐳)`$ that depends on $`𝐳`$ only via the corresponding value of the partition residue. It satisfies the following equations:

$$\left[1\tau G_0(\lambda )c(\mathrm{\Omega })\right]v(\mathrm{\Omega })=\frac{\tau \mathrm{\Phi }\mathrm{\hspace{0.17em}2}^n}{\lambda \alpha V_0}+\tau _{\mathrm{}}^{\mathrm{}}𝑑\mathrm{\Omega }^{}v(\mathrm{\Omega }^{})c(\mathrm{\Omega }^{})\chi (\mathrm{\Omega }^{},\mathrm{\Omega },\lambda ),$$
(43)
$$\chi (\mathrm{\Omega }^{},\mathrm{\Omega },\lambda )=\underset{r=1}{\overset{n}{}}\left(\genfrac{}{}{0pt}{}{n}{r}\right)G_r(\lambda )P_r(\mathrm{\Omega }^{}|\mathrm{\Omega }).$$
(44)
where $`\mathrm{\Phi }`$ is given in (35) and function $`c(x)`$ takes a set of discrete values (26). Using (35),(42) and (43) we obtain equations for $`u_𝐳`$

$$\left[1\tau G_0(\lambda )\epsilon _k\right]u_𝐳=\tau \underset{k^{}=1}{\overset{M}{}}\epsilon _k^{}\underset{𝐳^{}_k}{}G_{D(𝐳,𝐳^{})}(\lambda )u_𝐳^{}+\tau \epsilon _0\underset{𝐰_0}{}G_{D(𝐳,𝐰)}(\lambda )\varphi _𝐰,𝐳_k.$$
(45)
Decomposition (42) is only applied to amplitudes $`\varphi _𝐳`$ with $`𝐳_0`$. The system of equations for the components $`v(\mathrm{\Omega })`$ and $`u_𝐳`$ is closed by adding Eq. (35) for the amplitudes $`\varphi _𝐰`$ with $`𝐰_0`$ (ground states of the final Hamiltonian $`H_P`$) and taking (42) into account. We note that Eq.(43) for $`v(\mathrm{\Omega })`$ is coupled to the rest of the equations only via the symmetric term $`\mathrm{\Phi }`$

$`\mathrm{\Phi }=\overline{\mathrm{\Phi }}+\stackrel{~}{\mathrm{\Phi }}+\mathrm{\Phi }_0`$ (46)
$`\overline{\mathrm{\Phi }}=2^n{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xP(x)v(x)c(x),`$ (47)
$`\stackrel{~}{\mathrm{\Phi }}={\displaystyle \underset{k=1}{\overset{M}{}}}\epsilon _k{\displaystyle \underset{𝐳_k}{}}u_𝐳,\mathrm{\Phi }_0=\epsilon _0{\displaystyle \underset{𝐰_0}{}}\varphi _𝐰,`$
where distribution $`P(\mathrm{\Omega })`$ is given in (6).

#### IV.2.1  Minimum gap estimate for $`\omega _ME`$

We will analyze the above system of equations (42)-(47) assuming that the cutoff frequency $`\omega _M`$ satisfies Eq.(29). This condition corresponds to the linear region in the plot of the cumulative density of states given in insert to the Fig. II.1. According to Eqs. (6),(19), in this range the distribution functions $`P(\mathrm{\Omega })\overline{P}_{n/2}`$ and $`P_r(\mathrm{\Omega }^{}|\mathrm{\Omega })\overline{P}_r`$ take nearly constant values and spectral function $`\chi (\mathrm{\Omega }^{},\mathrm{\Omega },\lambda )`$ equals

$$\chi (\mathrm{\Omega }^{},\mathrm{\Omega },\lambda )\frac{s(\lambda )}{\sqrt{2\pi n\sigma ^2(0)}}.$$
(48)
where $`s(\lambda )`$ is given in (37). In this approximation, we can compute $`\stackrel{~}{\mathrm{\Phi }}`$ using equations for $`u_𝐳`$ in (45) and also the relations in (36), (37)

$$\stackrel{~}{\mathrm{\Phi }}=\kappa \left(\tau \mu s(\lambda )\right)\mathrm{\Phi }_0,\kappa (x)=\frac{x}{1+x}.$$
(49)
In the initial stage of the algorithm the amplitudes $`\varphi _𝐰`$ of the “solution” states are small $`|\mathrm{\Phi }_0|=𝒪(2^{n/2})`$. According to (49), we also have $`|\stackrel{~}{\mathrm{\Phi }}|=𝒪(2^{n/2})`$. Neglecting these terms and setting $`\mathrm{\Phi }\overline{\mathrm{\Phi }}`$, Eq. (43) gives a closed-form algebraic equation for $`\lambda `$

$$1+2\tau \mu \left(\frac{1}{\lambda \alpha V_0}+s(\lambda )\right)=0.$$
(50)
Expanding in a small parameter $`\mu 1`$ (cf.(29),(38)), we obtain the eigenvalue

$$\lambda _0^i(\tau )\alpha (\tau )V_02\tau \mu \frac{2(\tau \mu )^2\mathrm{ln}2}{\alpha }+𝒪(\mu ^3)(\alpha \mu ),$$
(51)
that accurately tracks the adiabatic ground state energy, $`\lambda _0(\tau )`$, from $`\tau =0`$, up until small vicinity of the avoided-crossing, $`\tau \tau ^{}`$ (see below) where $`|\mathrm{\Phi }_0|1`$.

In the avoided-crossing region, branch $`\lambda _0^i(\tau )`$ intersects with another branch, $`\lambda _0^f(\tau )`$, that tracks $`\lambda _0(\tau )`$ in the interval of time following the avoided-crossing, $`\tau ^{}<\tau 1`$. This branch corresponds to $`\overline{\mathrm{\Phi }}\mathrm{\Phi }_0,\stackrel{~}{\mathrm{\Phi }}`$. It can be obtained from simultaneous solution of equations for $`u_𝐳`$ (45) and $`\varphi _𝐰`$ that are approximately decoupled from Eq. (43) after $`\overline{\mathrm{\Phi }}`$ is neglected. Keeping this term in (45) gives rise to repulsion between branches $`\lambda _0^{i,f}(\tau )`$ at $`\tau =\tau ^{}`$ that determines the minimum gap $`g_{\mathrm{min}}`$ (see below).

To proceed, we obtain the equation for $`\mathrm{\Phi }_0`$ by adding equations for amplitudes $`\varphi _𝐰`$ that correspond to different states $`𝐰_0`$ and neglecting the coupling between these states separated by large Hamming distances, $`D(𝐰,𝐰^{})n/2`$. It can be shown using Eqs. (35) and (41)-(45) that $`u_𝐳`$ enters equation for $`\mathrm{\Phi }_0`$ through the term

$$\tau ^2\epsilon _0\underset{𝐳_0}{}_𝐳f_{𝐳,0}(\lambda )u_𝐳,$$
(52)
which is is a self-energy term corresponding to elementary bit-flip processes with initial and final states belonging to the subspace $`_0`$ (loop diagrams).

To express $`u_𝐳`$ in (52) through $`\varphi _𝐰`$ we solve Eq. (45) using order-by-order expansion in a small parameter $`n^1`$ (cf. Eqs. (36)-(41) and discussion there). In particular, one can show that to the leading order in $`n^1`$ the self-energy term (52) is determined by lowest-order loops with two bit flips that begin and end at $`_0`$. Then after some transformations, the equation for $`\mathrm{\Phi }_0`$ takes the form

$$\mathrm{\Phi }_0\left(\lambda \tau \epsilon _0\frac{\tau \alpha ^2\epsilon _0}{\lambda }\underset{𝐳^{}_0}{}\frac{\delta _{1,D(𝐳^{},𝐰)}}{\lambda \tau _𝐳^{}}\right)=\lambda \epsilon _0\tau d_02^n\left(\frac{\mathrm{\Phi }}{\lambda \alpha V_0}+\overline{\mathrm{\Phi }}s(\lambda )\right).$$
(53)
Here $`\alpha =1\tau `$ (cf. (34) and $`\overline{\mathrm{\Phi }}`$ is defined above. We now solve Eq. (53) jointly with (43) and obtain a closed-form equation for $`\lambda `$. We give it below in the region of interest $`|\tau 1/2|1`$

$`\left(\lambda \lambda _0^\mathrm{i}(\tau )\right)\left(\lambda \lambda _0^\mathrm{f}(\tau )\right)=n^22^n\mathrm{\Delta }^2/4`$ (54)
$`\mathrm{\Delta }d_0^{1/2}\left(1+\mu \tau ^{}\mathrm{ln}2+𝒪(\mu ^2)\right),`$
where the branch $`\lambda _0^\mathrm{i}(\tau )`$ is given above and the branch $`\lambda _0^\mathrm{f}(\tau )`$ satisfies Eq. (53) with r.h.s. there set to zero,

$$\lambda _0^\mathrm{f}(\tau )\tau \epsilon _01/2,|\tau 1/2|1.$$
(55)
Avoided-crossing in (54) takes place at $`\tau =\tau ^{}`$

$$\lambda _0^\mathrm{i}(\tau ^{})=\lambda _0^\mathrm{f}(\tau ^{}),\tau ^{}\frac{1}{2}+\frac{1}{4n}\mathrm{log}_2\frac{d_0}{\mu }.$$
(56)
The value of minimum gap between the two roots of (54) equals

$$g_{\mathrm{min}}=n\mathrm{\Delta }\mathrm{\hspace{0.17em}2}^{n/2}.$$
(57)
where $`\mathrm{\Delta }`$ is defined in (27).

Based on the above analysis one can also estimate the matrix element $`|\varphi _1|\stackrel{~}{H}_\tau |\varphi _0|_{\tau =\tau ^{}}n`$. Then from Eq. (25) (see also discussion after Eq. (28)) one can estimate the run-time of the quantum adiabatic algorithm

$$T\frac{d_0|H_{\tau 01}^{}|}{g_{\mathrm{min}}^2}=𝒪((nd_0)^12^n).$$
(58)
It follows from the above that eigenvalue branch $`\lambda _0^\mathrm{i}(\tau )`$ corresponds to a state,

$$|\varphi _0\underset{𝐳\{0,1\}^n}{}v(\mathrm{\Omega }_z)|𝐳,$$

which is extended in the space of the bit configurations $`|𝐳`$: according to (43) it contains a large number ($`𝒪(2^n)`$) of exponentially small ($`𝒪(2^{n/2})`$) individual amplitudes. This state originates at $`\tau =0`$ from the totally symmetric initial state $`|\psi (0)`$ (23). In the small region $`|\tau \tau ^{}|g_{\mathrm{min}}`$ it is transformed into the state that corresponds to the eigenvalue branch $`\lambda _0^\mathrm{f}(\tau )`$ and is localized in Hamming distances $`D(𝐳,𝐰)`$ near the subspace $`𝐰_0`$ containing the solution to NPP $`|\varphi _0_{𝐰_0}|𝐰`$. Minimum gap at the avoided-crossing is determined by the overlap between the extended and localized states.

At later times $`\tau >\tau ^{}`$ a similar picture applies to the avoided crossing of the extended-state energy $`\lambda _0^\mathrm{i}(\tau )`$ with energies of localized states $`\lambda _k^\mathrm{f}(\tau )`$ corresponding to $`𝐳_k`$ with $`1kn`$ (excited levels of the final Hamiltonian $`H_P`$ (28)). The existence of the extended eigenstate of $`\stackrel{~}{H}(\tau )`$ whose properties do not depend on a particular instance of NPP follows directly from Eq. (43) that involves only a self-averaging quantity $`\chi (\mathrm{\Omega }^{},\mathrm{\Omega },\lambda )`$. This quantity varies smoothly over the broad range of partition residues $`|\mathrm{\Omega }^{}|,|\mathrm{\Omega }|E`$ and does not allow for the compression of the wave-packet $`v(\mathrm{\Omega }_z)`$ on the much smaller scale $`𝒪(2^n)`$. This gives rise to an eigenstate with probability amplitude of individual states $`|𝐳`$ that depends smoothly on energy in this range.

#### IV.2.2 Analysis of the general case

The above picture of avoided-crossing remains qualitatively the same when the condition (29) is relaxed (cf. insert in the Fig. 2). Away from the avoided-crossing point, $`\tau <\tau ^{}`$, the ground state wavefunction $`v(\mathrm{\Omega }_𝐳)`$ and energy $`\lambda _0^i(\tau )`$ are obtained directly from Eq. (43) with replacement $`\mathrm{\Phi }\overline{\mathrm{\Phi }}`$ and Eq. (44) taken into account. Because the spectral function $`\chi (\mathrm{\Omega },\mathrm{\Omega }^{},\lambda )`$ changes only slightly on the scale $`E_{\mathrm{min}}=𝒪(n^{1/2}2^n)`$ the wave packet $`_𝐳v(\mathrm{\Omega }_𝐳)|𝐳`$ remains extended, $`|v(\mathrm{\Omega }_𝐳|=𝒪(2^{n/2})`$, and therefore $`\mathrm{\Phi }_0=𝒪(2^{n/2})`$.

Beyond the avoided-crossing point, $`\tau >\tau ^{}`$, the ground state is localized near $`𝐰`$ and eigenvalue branch $`\lambda _0^f(\tau )`$ is obtained from Eq. (53) with r.h.s. set to zero (cf. Sec. IV.2.1). The point $`\tau =\tau ^{}`$ is located at the intersection of the two branches $`\lambda _0^i(\tau )\lambda _0^f(\tau )`$ and the level repulsion is of the order of the overlap factor between the extended and localized states

$$g_{\mathrm{min}}\underset{𝐰_0}{}v(\mathrm{\Omega }_𝐰)2^{n/2}.$$
(59)
Ground-state wavefunction $`\varphi _𝐳`$ at the avoided-crossing is shown in Fig. 2 for modest value of $`n`$, but the separation into slowly- and rapidly-varying parts (42) is clearly seen.

We did not perform a direct numerical study of the dependence of $`g_{\mathrm{min}}`$ on $`n`$ since we only simulated adiabatic eigenvalues for small instances of NPP. We argue, however, that even for a fixed $`n`$ the scaling of $`g_{\mathrm{min}}`$ with $`n`$ can be inferred from the shape of the cumulative density of states

$$\eta (\lambda )=_0^\lambda 𝑑x\underset{k=0}{\overset{k_m}{}}\delta \left(\lambda _kx\right),k_m=2^n1,$$
(60)
where $`\lambda _k\lambda _k(\tau )`$ are eigenvalues of $`H(\tau )`$ (24). These eigenvalues are plotted in Fig. 3 near the avoided-crossing $`\tau =\tau ^{}`$ where the spectrum of $`\lambda _k`$ is quasi-continuous. The shape of the plot is well approximated by the square-root function:

$$\lambda _\eta \mathrm{const}+\left(\frac{\eta }{\eta _m}\right)^{1/2},\eta _m=𝒪(2^n).$$
(61)
It is clear that for $`\eta 1`$ we have $`\lambda _\eta 2^{n/2}`$ which corresponds to Eq. (57). Note that this qualitative analysis is based on the assumption that the asymptotic properties of $`\lambda _0`$ for large $`n`$ can be inferred from the behavior of $`\lambda _\eta `$ for $`\eta 1`$.

## V  Simulations of time-dependent Schrödinger equation

We also study the complexity of the algorithm by numerical integration of the time-dependent Schrödinger equation with Hamiltonian $`H(t)`$ and initial state $`|\psi (0)`$ defined in Eqs. (22),(23),(27),(28). Here we relax the condition $`\omega _ME`$ used above in the analytical treatment of the problem; in simulations the value of $`M`$ is set automatically to be an integer closest to $`\mathrm{log}_2_{j=0}^na_j`$ (cf. (27)). We introduce a complexity metric for the algorithm, $`C(T)=(1+T)d_0/p_0(T)`$ where $`p_0(t)=_{𝐰_0}|\psi _𝐰(t)|^2`$. A typical plot of $`C(T)`$ for an instance of the problem with $`n`$=15 numbers is shown in the insert of Fig. 4. At very small $`T`$ the wavefunction is close to the symmetric initial state and the complexity is $`2^n`$. The extremely sharp decrease in $`C(T)`$ with $`T`$ is due to the buildup of the population $`p_0(T)`$ in the ground level, $`_𝐳=\epsilon _0`$, as quantum evolution approaches the adiabatic limit. At certain $`T=T_{\mathrm{min}}`$ the function $`C(T)`$ goes through the minimum: for $`T>T_{\mathrm{min}}`$ the decrease in the number of trials $`d_0/p_0(T)`$ does not compensate anymore for the overall increase in the runtime $`T`$ for each trial. For a given problem instance the “minimum” complexity $`C_{\mathrm{min}}=C(T_{\mathrm{min}})`$ is obtained via one dimensional minimization over $`T`$. The plot of the complexity $`C_{\mathrm{min}}`$ for different values of $`n`$ in Fig. 1 appears to indicate the exponential scaling law, $`C_{\mathrm{min}}2^{0.8n}`$ for not too small values of $`n`$ 11.

## VI Discussion

In conclusion, we have developed a general method for the analysis of avoided-crossing phenomenon in quantum spin-glass problems and used it to study the performance of the quantum adiabatic evolution algorithm on random instances of the Number Partitioning problem. This algorithm is viewed as a “quantum local search” with matrix elements of the Green function $`G_r`$ ($`r=1,\mathrm{},n1`$) giving the quantum amplitudes of the transitions with different number of spin flips $`r`$. Our approach is similar to the analysis of a quantum diffusion in a disordered medium with the model of disorder defined by the one- and two-point distribution functions $`P(\mathrm{\Omega }),P_{r,𝐳}(\mathrm{\Omega }^{})`$.

We have shown that the conditional distribution of partition residues $`P_r(\mathrm{\Omega }^{}|\mathrm{\Omega })`$ in the neighborhood of a given string formed by all possible $`r`$-bit flips depends on the value of the partition residue for that string but not on the string itself. This is a specific property of the random Number Partitioning problem.

We used the above property to describe a quantum diffusion in the energy space (Eq. (43)). This reduction in the dimensionality leads to the formation of the eigenstate which is extended in the energy space. Near the avoided-crossing the adiabatic ground state changes from extended to mostly localized near the solution to the optimization problem. Because the extended and localized state amplitudes are nearly orthogonal to each other the repulsion between the corresponding branches of eigenvalues (the minimum gap) is exponentially small, $`g_{\mathrm{min}}n\mathrm{\hspace{0.17em}2}^{n/2}`$, and the run time of the algorithm scales exponentially with $`n`$. Analytical results are in qualitative agreement with numerical simulations of the time-dependent Schrödinger equation for small-to-moderate instances of the Number Partitioning problem ($`n17`$).

One can show that the effect of quantum diffusion in reduced-dimensional space that leads to the formation of the extended state can also occur in other random NP-complete problems SST . The method developed in this paper will be applied to study the performance of continuous-time quantum algorithms for different random combinatorial optimization problems. Also the present framework can be applied to the analysis of quantum annealing algorithms for combinatorial optimization problems kadowaki ; santoro . This is a classical algorithm that is conceptually very close to the quantum adiabatic evolution algorithm considered above kadowaki1 . The former uses the Quantum Monte Carlo method to simulate on classical computers a partition function and ground-state energy of a quantum system with slowly varying Hamiltonian that merges at the final moment with the problem Hamiltonian of a given classical optimization problem. Among other possible applications of our method is the analysis of tunneling phenomenon in the low-temperature dynamics of random magnets.

We note that the specific property of the Number Partitioning problem (that distinguishes it from the other NP-complete problems) is a very weak dependence of $`P_r(\mathrm{\Omega }^{}|\mathrm{\Omega })`$ on $`\mathrm{\Omega }`$ for not too large values of $`\mathrm{\Omega }^{},\mathrm{\Omega }\sqrt{r(nr)}`$ that takes place for all values of $`r[1,n1]`$. This rapid fall-off of correlations during the local search (both classical and quantum) is a reason that the exponential complexity of optimization algorithms for the Number Partitioning problem can be seen already for the relatively small values of $`n15`$ (cf. Fig.4).

Finally, our analysis of sub-harmonic resonances in the Fourier transform $`I(s)`$ of the distribution function $`P(\mathrm{\Omega })`$ suggests a possible connection between NPP and the integer factorization problem. If, for a given set of $`a_j`$’s, there is a number $`q`$ that satisfies the condition (64) then dividing all numbers $`a_j`$ by $`q`$ we obtain a new instance of NPP with numbers $`k_j=a_j/q`$ that will be completely equivalent to the old one. It is important that the precision of the numbers $`k_j`$ is restricted by $`b\mathrm{log}_2q`$. If the value of $`q`$ is sufficiently large, $`\mathrm{log}_2qbn`$, then $`k_j`$’s correspond to a low precision instance of NPP, i.e. to the computationally easy phase mentioned in Sec. II. This is exactly the case when sub-harmonic resonances become substantial. One can fix the parameter $`\xi =b/n1`$ in a high-precision (computationally hard) case and compute, for randomly generated instances $`\{a_j\}`$ an approximate greatest common divider, i.e. a largest number $`q`$ that satisfies (64). The distribution of these numbers determines a fraction of high-precision instances of NPP (out of all possible $`2^{nb}`$ problem instances) that really belong to a low-precision (computationally easy) “phase”.

Advance knowledge of this information would be of importance if one is using NPP for encryption purposes cr , especially because NPP is otherwise a very difficult problem for both quantum and classical computers Mertens2000 . It is not obvious at this stage what the asymptotic form of this distribution will be in the limit of large $`n`$ (cf. Fig. A).

We are not aware of any classical algorithm that could verify if such a number $`q`$ exists for a given set of $`a_j`$ in a time polynomial in both $`n`$ and $`b`$. However, on a quantum computer one can apply a Shor algorithm to test in polynomial time if strong sub-harmonic resonances exist. This question is deferred to a future study.

## VII Acknowledgments

The authors benefited from stimulating discussions with P. Cheeseman, R.D. Morris (NASA ARC) and U. Vazirani (UC Berkley). We also acknowledge the help of J. Lohn (Automated Design of Complex Systems group, NASA ARC) for providing computer facilities. This research was supported by NASA Intelligent Systems Revolutionary Computing Algorithms program (project No: 749-40), and also by NASA Ames NAS Center.

## Appendix A  Sub-harmonic resonances

We note that function $`I(s)`$ in (5) can also have additional sharp resonances in the range $`0<|s|2^b`$. To understand their origin we consider first a particular case when rational $`b`$-bit numbers $`a_1,a_2,\mathrm{},a_n`$ all have a number $`q>`$ 2<sup>-b</sup> as a “common divisor”, i.e., there exist integers $`k_1,k_2,\mathrm{},k_n`$ such that

$$\frac{a_1}{k_1}=\frac{a_2}{k_2}=\mathrm{}=\frac{a_n}{k_n}=q.$$
(62)
In this case additional resonances of $`I(s)`$ occur at the multiples of $`\pi /q`$. Assume now that $`q`$ is no longer an exact divisor of numbers $`a_j`$ but all the residues of the divisions $`a_j/q`$ are sufficiently small. Then contributions from the additional resonances at $`sm\pi /q`$ ($`m=1,\mathrm{\hspace{0.17em}2},\mathrm{}`$) to the integral in (5) can be estimated as follows (for simplicity we give a result for the case $`En^{1/2}`$):

$`P(0){\displaystyle \frac{2^n}{\sqrt{2\pi n\sigma ^2(0)}}}e^{\gamma (q)}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\zeta \left({\displaystyle \frac{m\pi \eta }{2q}}\right)(1)^{mp}`$ (63)
$`p={\displaystyle \underset{j=1}{\overset{n}{}}}\left[{\displaystyle \frac{a_j}{q}}\right],\gamma (q)={\displaystyle \frac{\pi ^2}{2}}\left(\left\{{\displaystyle \frac{a_j}{q}}\right\}^2{\displaystyle \frac{a_j}{\sqrt{\pi n\sigma ^2(0)}}}\left\{{\displaystyle \frac{a_j}{q}}\right\}\right)`$
Here $`[x]`$ and $`\{x\}`$ denote integer and fractional parts of a number $`x`$, respectively. If the total “dephasing” factor $`e^{\gamma (q)}1`$, then contribution (63) cannot be neglected in the steepest-descent analysis of (5) (in general, on should keep contributions from all divisors $`q`$ with small dephasing factors $`e^{\gamma (q)}`$).

We note that the window function $`\zeta \left(\frac{m\pi \eta }{2q}\right)1`$ for $`q2^n`$ and it decays to zero at smaller values of $`q`$. We studied numerically the greatest root $`q_{\mathrm{max}}`$ of the the algebraic equation

$$\gamma (q)=\gamma _c$$
(64)
for a fixed value of $`\gamma _c1`$. For the sets of random $`b`$-bit numbers $`a_j`$ the dependence of the mean value of $`q_{\mathrm{min}}`$ on the problem size $`n<b`$ is shown in Fig. 5. For $`nb`$ we have exponential decrease of $`q_{\mathrm{max}}`$ with $`n`$ and for larger values of $`nb`$ the value of $`q_{\mathrm{min}}`$ steeply drops to 1. According to the discussion above, in order to neglect the saddle-points with $`s>0`$ in (5) (additional resonances) the value of $`q_{\mathrm{min}}`$ should satisfy the following condition in the asymptotic limit $`b\mathrm{}`$:

$$q_{\mathrm{max}}\mathrm{max}[2^n,2^b],1nb,$$
(65)
with $`\gamma _c`$ fixed at some small constant value. Because the precision $`b`$ that we used in the simulations was not very high (limited by machine precision) it is not possible to obtain the asymptotic form of the dependence of $`q_{\mathrm{min}}`$ on $`n`$ in the range given in (65). Neither we can describe the shape of the plot in Fig. 5 analytically in that range. However, it appears from the figure that the condition (65) is satisfied for sufficiently large $`n`$.

## Appendix B Properties of the conditional distribution of signed residues in NPP

We perform the summation over the spin configurations in Eq. (17) with Eq. (15) taken into account. Similar to the derivation of Eq. (5) we use integral representation for delta function and obtain

$$P_r(\mathrm{\Omega },\mathrm{\Omega }^{})=\left(\genfrac{}{}{0pt}{}{n}{r}\right)^1_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\frac{dsds^{}}{4\pi ^2}\zeta \left(\frac{\mathrm{\Delta }\mathrm{\Omega }s}{2}\right)\zeta \left(\frac{\mathrm{\Delta }\mathrm{\Omega }^{}s^{}}{2}\right)\underset{𝐉}{}U_𝐉(s,s^{}),$$
(66)
$$U_𝐉(s,s^{})=\underset{j𝐉}{}\mathrm{cos}(a_j(ss^{}))\times \underset{j𝐉}{}\mathrm{cos}(a_j(s+s^{}))\times e^{i(s\mathrm{\Omega }+s^{}\mathrm{\Omega }^{})}.$$

Here the sum is over all possible subsets $`𝐉=\{j_1,j_2,\mathrm{},j_r\}`$ of length $`r`$ obtained from the set of integers $`j=1,2,\mathrm{},n`$. Window function $`\zeta (x)`$ is defined in (5). After the change of variables

$$x^{}=s+s^{},x=ss^{},$$
(67)
we obtain from (66) that $`U_𝐉(s,s^{})`$ factorizes into a product of two terms

$$U_𝐉(s,s^{})=𝒱_𝐉(x)\stackrel{~}{𝒱}_𝐉(x^{})$$

$$𝒱_𝐉(x)=\mathrm{exp}\left(\frac{ix(\mathrm{\Omega }\mathrm{\Omega }^{})}{2}\right)\underset{j𝐉}{}\mathrm{cos}(a_jx),\stackrel{~}{𝒱}_𝐉(x^{})=\mathrm{exp}\left(\frac{ix^{}(\mathrm{\Omega }+\mathrm{\Omega }^{})}{2}\right)\underset{j𝐉}{}\mathrm{cos}(a_jx^{}).$$
(68)
In what follows we will analyze several limiting cases.
$`r,nr1`$:
In this case both functions $`𝒱_𝐉(x)`$ and $`\stackrel{~}{𝒱}_𝐉(x^{})`$ are very steep and similar to the analysis in Sec.II.1 integrals in (66) can be evaluated by the steepest descent method. With the appropriate choice of the coarse-graining windows $`\mathrm{\Delta }\mathrm{\Omega },\mathrm{\Delta }\mathrm{\Omega }^{}`$ in (66) (see below) contribution to the integrals comes from the vicinity of the point ($`x=0,x^{}=0`$). Near this point we use

$$\underset{j𝐉}{}\mathrm{cos}(a_jx)\mathrm{exp}\left(\frac{r(x\sigma _𝐉)^2}{2}\right),\underset{j𝐉}{}\mathrm{cos}(a_jx)\mathrm{exp}\left(\frac{(nr)(x^{}\overline{\sigma }_𝐉)^2}{2}\right)$$
(69)
where

$$(\sigma _𝐉)^2=\frac{1}{r}\underset{j𝐉}{}a_j^2,(\overline{\sigma }_𝐉)^2=\frac{1}{nr}\underset{j𝐉}{}a_j^2.$$

Since each sum here contains a large number of terms we obtain for i.i.d. random numbers $`a_1,\mathrm{},a_n`$ (cf. (6))

$$(\sigma _𝐉)^2\sigma ^2(0)+𝒪\left(\frac{1}{r}\right),(\overline{\sigma }_𝐉)^2\sigma ^2(0)+𝒪\left(\frac{1}{nr}\right),$$
(70)
where $`\sigma ^2=a^2`$ is given in (6). Using Eqs. (68)-(70) and replacing the window functions in (66) by unity, we compute the Gaussian integrals in (66) and obtain

$$P_r(\mathrm{\Omega },\mathrm{\Omega }^{})=\frac{1}{4\pi \sigma ^2(0)\sqrt{r(nr)}}\mathrm{exp}\left[\frac{1}{8\sigma ^2(0)}\left(\frac{(\mathrm{\Omega }\mathrm{\Omega }^{})^2}{r}+\frac{(\mathrm{\Omega }+\mathrm{\Omega }^{})^2}{nr}\right)\right].$$
(71)
The size of the coarse-graining windows in (66) is chosen to satisfy the conditions

$$2^n\left(\genfrac{}{}{0pt}{}{n}{r}\right)^1\mathrm{\Delta }\mathrm{\Omega }\mathrm{\Delta }\mathrm{\Omega }^{}\sqrt{r(nr)}$$

From Eq. (71) and Eq. (6) one can directly obtain the conditional distribution function $`P_r(\mathrm{\Omega }|\mathrm{\Omega }^{})`$

$$P_r(\mathrm{\Omega }^{}|\mathrm{\Omega })=\frac{1}{\sqrt{2\pi n\sigma ^2(q)}}\mathrm{exp}\left[\frac{\left(\mathrm{\Omega }^{}q\mathrm{\Omega }\right)^2}{2n\sigma ^2(q)}\right].$$
(72)
$`r=1`$; $`r=n1`$:
For $`r=1`$ function $`\overline{𝒱}_𝐉(x^{})`$ contains a product of $`n1`$ terms and is very steep. The corresponding integral over $`x^{}`$ in (66) should be taken by the steepest descent method. However $`𝒱_𝐉(x)`$ simply oscillates at frequencies $`(\mathrm{\Omega }\mathrm{\Omega }^{})/2\pm a_j`$ and the integral over $`x`$ in (66) should be evaluated using the corresponding oscillating factors. In the opposite case $`r=n1`$, function $`𝒱_𝐉(x)`$ is very steep and the integral over $`x`$ in (66) should be taken by steepest descent. But the integral over $`x^{}`$ there should be evaluated using $`\overline{𝒱}_𝐉(x^{})`$ that oscillates at the frequencies, $`(\mathrm{\Omega }+\mathrm{\Omega }^{})/2\pm a_j`$. Finally, one can obtain using i.i.d. numbers $`a_j`$’s in $`[0,1]`$ interval :

$$P_r(\mathrm{\Omega }^{}|\mathrm{\Omega })=\frac{1}{4}\left[\mathrm{\Theta }\left(\mathrm{\Omega }\mathrm{\Omega }^{}+2\right)\mathrm{\Theta }\left(\mathrm{\Omega }\mathrm{\Omega }^{}2\right)\right]+𝒪\left(\frac{1}{n}\right),(r=1,n1).$$
(73)
The minus (plus) sign in (73) corresponds to $`r=1`$ ($`r=n1`$). Similarly one can obtain the result for any fixed value of $`r`$ or $`nr`$ (that does not scale with $`n`$). For $`|\mathrm{\Omega }|,|\mathrm{\Omega }^{}|1`$ (73) is reduced to (19).

### Numerical simulations of conditional distribution $`P_{r,𝐳}(\mathrm{\Omega }^{})`$

We compute the following integrated quantity:

$$Q=\frac{1}{2}_0^\mathrm{\Omega }^{}𝑑\eta P_{r,𝐳}(\eta ),$$
(74)
for different values of $`r`$, $`\mathrm{\Omega }^{}`$ and different strings $`𝐳`$ with $`E_𝐳1`$. Numerical results are compared in the insert to Fig.1 with theoretical result below obtained using $`P_r(\mathrm{\Omega }^{}|\mathrm{\Omega })`$ from Eq. (72)

$$\frac{1}{2}_0^\mathrm{\Omega }^{}𝑑\eta P_r(\eta |0)=\mathrm{erf}\left(\frac{\mathrm{\Omega }^{}}{\sigma (q)\sqrt{2n}}\right).$$
(75)
Theoretical and numerical curves nearly coincide with each other. To accurately compare the normalization factor in (72) (see also (20)) we compare the theoretical results with numerical values of $`P_{r,𝐳}(0)`$ for different $`r`$ and strings $`𝐳`$ corresponding to $`E_𝐳1`$. The results are plotted in Fig. 1.


format: text

text:
# On the Multiboson Coherent State in Deformed Boson Scheme

Recently, the present authors published four papers concerning the deformed boson scheme including conventional $`q`$-deformation in time-dependent variational method. The series of these papers consists of (I), (II), (III) and the short note and, in this series, we discussed various aspects of the $`q`$-deformation appearing in the boson coherent state and its generalized form. For example, in (I), a possible form of the multiboson coherent state and its generalization were described. Its treatment is, in some sense, in close connection with that given by Penson and Solomon. However, the treatment in (I) is limited to a certain special boson subspace. In this sense, it is unsatisfactory and to make it complete is the aim of this paper. Further, the present paper gives a generalization of the treatment proposed by the present authors, which is related to a new boson realization in the Lipkin model.

In (I), a possible form of the multiboson coherent state, which we denote as $`|c_0`$, was defined in the boson space $`\stackrel{~}{B}`$ composed by the boson operator $`(\stackrel{~}{c},\stackrel{~}{c}^{})`$ in the following form :

$$|c_0=\left(\sqrt{\mathrm{\Gamma }_0}\right)^1\underset{n=0}{\overset{\mathrm{}}{}}\left(\sqrt{(mn)!}\right)^1\gamma ^n|mn.$$
(1)
Here, $`\gamma `$ denotes a complex parameter. The quantity $`\mathrm{\Gamma }_0`$ is the normalization constant given as

$$\mathrm{\Gamma }_0=\underset{n=0}{\overset{\mathrm{}}{}}\left[(mn)!\right]^1(|\gamma |^2)^n.$$
(2)
The state $`|mn`$ is expressed in terms of

$`|mn=\left(\sqrt{(mn)!}\right)^1(\stackrel{~}{c}^{})^{mn}|0,(n=0,1,2,\mathrm{})`$ (3a)
$`\stackrel{~}{c}|0=0.`$ (3b)
The operator $`(\stackrel{~}{c}^{})^m`$ is the building block of the state $`|c_0`$ and, in this sense, we can call $`|c_0`$ the multiboson coherent state in which $`m`$ denotes an integer ($`m=2,3,\mathrm{}`$) characterizing the multiboson. Clearly, the set $`\{|mn\}`$ is an orthogonal but not a complete one : The set composes a subspace of $`\stackrel{~}{B}`$, which we call $`\stackrel{~}{B}_0`$.

With the aim of investigating the structure of the deformation of the multiboson coherent state, in (I), we adopted the basic idea of the MYT boson mapping. For the image of the space $`\stackrel{~}{B}_0`$, we prepare a boson space $`\widehat{B}_0`$ constructed by the boson operator $`(\widehat{c},\widehat{c}^{})`$. The orthogonal set of $`\widehat{B}_0`$ consists of the states

$`|n=(\sqrt{n!})^1(\widehat{c}^{})^n|0,(n=0,1,2,\mathrm{})`$ (4a)
$`\widehat{c}|0=0.`$ (4b)
The set $`\{|n\}`$ is complete. Following the basic idea of the MYT boson mapping, we make the set $`\{|mn\}`$ correspond to $`\{|n\}`$ :

$$|mn|n,\text{i.e.,}|n=𝑼_0|mn.$$
(5)
The operator $`𝑼_0`$ plays a role of the mapping operator from $`\stackrel{~}{B}_0`$, a subspace of $`\stackrel{~}{B}`$, to $`\widehat{B}_0`$ and it is given by

$$𝑼_0=\underset{n=0}{\overset{\mathrm{}}{}}|nmn|.$$
(6)
The operator $`𝑼_0`$ satisfies the relation

$$𝑼_0𝑼_0^{}=1,𝑼_0^{}𝑼_0=\stackrel{~}{P}_m.(\stackrel{~}{P}_m=\underset{n=0}{\overset{\mathrm{}}{}}|mnmn|)$$
(7)
Then, the image of $`|c_0`$ mapped from $`\stackrel{~}{B}_0`$ to $`\widehat{B}_0`$, which we denote $`|c_0`$, is given in the following form :

$`|c_0`$ $`=`$ $`𝑼_0|c_0`$ (8)
$`=`$ $`\left(\sqrt{\mathrm{\Gamma }_0}\right)^1{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}f_m^0(n)(\sqrt{n!})^1\gamma ^n|n,`$
$`f_m^0(n)`$ $`=`$ $`\sqrt{n!}\left(\sqrt{(mn)!}\right)^1.`$ (9)
The form (8) shows that $`|c_0`$ is deformed from the conventional boson coherent state by the function $`f_m^0(n)`$. This means that by changing $`f_m^0(n)`$ from the form (9), we have various deformed states. The state $`|c_0`$ can be rewritten as

$`|c_0=\left(\sqrt{\mathrm{\Gamma }_0}\right)^1\mathrm{exp}\left((\sqrt{m^m})^1\gamma \widehat{c}^{}\left(\sqrt{F_m^0(\widehat{N}_c)}\right)^1\right)|0,`$ (10)
$`F_m^0(\widehat{N}_c)={\displaystyle \underset{p=1}{\overset{m1}{}}}(\widehat{N}_c+p/m)=(\widehat{N}_c+1)^1{\displaystyle \underset{p=1}{\overset{m}{}}}(\widehat{N}_c+p/m).\widehat{N}_c=\widehat{c}^{}\widehat{c}.`$ (11)
The form (10) for $`m=2`$ can be seen in Appendix of Ref.?. The state $`|c_0`$ is an eigenstate of $`\widehat{\gamma }`$ defined in the following :

$$\widehat{\gamma }|c_0=\gamma |c_0,\widehat{\gamma }=\sqrt{m^m}\sqrt{F_m^0(\widehat{N}_c)}\widehat{c}.$$
(12)
The form (10) also shows that $`|c_0`$ is deformed from the conventional boson coherent state. However, the above treatment contains unsatisfactory point. The system is treated only in the subspace $`\stackrel{~}{B}_0`$ of $`\stackrel{~}{B}`$. Therefore, in order to complete the treatment, it may be necessary to investigate the effect coming from the subspace orthogonal to $`\stackrel{~}{B}_0`$. The aim of this note is to investigate the above-mentioned problem.

For our present aim, we treat the system in the whole space $`\stackrel{~}{B}`$. In this case, the orthogonal set consists of the following states :

$`|mn+r=\left(\sqrt{(mn+r)!}\right)^1(\stackrel{~}{c}^{})^{mn+r}|0.`$
$`(n=0,1,2,\mathrm{};r=0,1,2,\mathrm{},m1)`$ (13)
Of course, the set $`\{|mn+r\}`$ is complete and the space spanned by the states with $`r=0`$ corresponds to the subspace $`\stackrel{~}{B}_0`$. The case $`m=2`$ was investigated in Ref.? in a form different from the present one.

In order to obtain the image of the space $`\stackrel{~}{B}`$, let us prepare a boson space consisting of not only the boson $`(\widehat{c},\widehat{c}^{})`$ but also the boson $`(\widehat{d},\widehat{d}^{})`$, the orthogonal set of which is given by

$$|n,r=|n|r),$$
(14)
$`|r)=(\sqrt{r!})^1(\widehat{d}^{})^r|0),(r=0,1,2,\mathrm{})`$ (15a)
$`\widehat{d}|0)=0.`$ (15b)
We call the above space as $`\overline{B}`$. The state $`|n`$ is defined in the relation (4a). We make the set $`\{|mn+r\}`$ correspond to the set $`\{|n,r\}`$ :

$$|mn+r|n,r,\text{i.e.,}|n,r=𝑼|mn+r.$$
(16)
Since $`r=0,1,2,\mathrm{},m1`$ in the set $`\{|mn+r\}`$, $`r`$ in the state $`|n,r`$ shown in the relation (16) is also equal to $`0,1,2,\mathrm{},m1`$. This means that the set $`\{|n,r\}`$ as the image of $`\{|mn+r\}`$ composes a subspace of $`\overline{B}`$, which we denote $`\widehat{B}`$. The operator $`𝑼`$ is written as

$$𝑼=\underset{n=0}{\overset{\mathrm{}}{}}\underset{r=0}{\overset{m1}{}}|n,rmn+r|.$$
(17)
The operator $`𝑼`$ satisfies the relation

$$𝑼^{}𝑼=1,𝑼𝑼^{}=\widehat{P}_m.(\widehat{P}_m=\underset{r=0}{\overset{m1}{}}|r)(r|)$$
(18)
As a possible and natural extension of the multiboson coherent state $`|c_0`$, we define the following state in the space $`\stackrel{~}{B}`$ :

$`|c=\left(\sqrt{\mathrm{\Gamma }}\right)^1{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{r=0}{\overset{m1}{}}}\left(\sqrt{(mn+r)!}\right)^1\gamma ^n\delta ^r|mn+r,`$ (19)
$`\mathrm{\Gamma }={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{r=0}{\overset{m1}{}}}[(mn+r)!]^1(|\gamma |^2)^n(|\delta |^2)^r.`$ (20)
Here, $`\delta `$ denotes a complex parameter additional to $`\gamma `$. If restricted to $`r=0`$ in the summation in (19) and (20), $`|c`$ is reduced to $`|c_0`$. With the use of $`𝑼`$ given in the relation (17), the image of $`|c`$, which we denote $`|c`$, can be derived :

$`|c`$ $`=`$ $`𝑼|c`$ (22)
$`=`$ $`\left(\sqrt{\mathrm{\Gamma }}\right)^1{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{r=0}{\overset{m1}{}}}f_m(n,r)(\sqrt{n!})^1(\sqrt{r!})^1\gamma ^n\delta ^r|n,r,`$
$`f_m(n,r)=\sqrt{n!}\sqrt{r!}\left(\sqrt{(mn+r)!}\right)^1.`$
The form (22) shows that the state $`|c`$ is different from the conventional boson coherent state. The state $`|c`$ can be rewritten in the form

$`|c=\left(\sqrt{\mathrm{\Gamma }}\right)^1\mathrm{exp}\left(\left(\sqrt{m^m}\right)^1\gamma \widehat{c}^{}\left(\sqrt{F_m(\widehat{N}_c,\widehat{N}_d)}\right)^1\right)\mathrm{exp}\left(\delta \widehat{d}^{}\widehat{P}_{m1}\right)|0.`$
$`(|0|0)=|0)`$ (23)
$`F_m(\widehat{N}_c,\widehat{N}_d)=\left(\widehat{N}_c+1\right)^1{\displaystyle \underset{p=1}{\overset{m}{}}}\left[\widehat{N}_c+(\widehat{N}_d+p)/m\right],\widehat{N}_d=\widehat{d}^{}\widehat{d},`$ (24)
$`\widehat{P}_{m1}={\displaystyle \underset{r=0}{\overset{m2}{}}}|r)(r|.`$ (25)
The function $`F_m(\widehat{N}_c,\widehat{N}_d)`$ is a natural extension from $`F_m^0(\widehat{N}_c)`$ shown in the relation (11). The part $`\mathrm{exp}(\delta \widehat{d}^{}\widehat{P}_{m1})|0`$ is a concrete example of the case investigated in Ref.?.

The state $`|c`$ has a close connection with the following operators :

$`\widehat{\gamma }=\sqrt{m^m}\sqrt{F_m(\widehat{N}_c,\widehat{N}_d)}\widehat{c},`$ (26)
$`\widehat{\delta }=\sqrt{(m\widehat{N}_c+\widehat{N}_d+1)(\widehat{N}_d+1)^1}\widehat{d}.`$ (27)
The connection can be seen in the relation

$`\widehat{\gamma }|c=\gamma |c,`$ (28)
$`\widehat{\delta }|c=\delta \widehat{P}_{m1}|c.`$ (29)
It is useful to note that the following relation for $`F_m`$ is satisfied :

$$F_m(\widehat{N}_c,\widehat{N}_d+1)=\left[\widehat{N}_c+1+(\widehat{N}_d+1)/m\right]\left[\widehat{N}_c+(\widehat{N}_d+1)/m\right]^1F_m(\widehat{N}_c,\widehat{N}_d).$$
(30)
With use of the relation (30), we can prove

$$[\widehat{\delta },\widehat{c}^{}\left(\sqrt{F_m(\widehat{N}_c,\widehat{N}_d)}\right)^1]=0.$$
(31)
Further, we have

$$[\widehat{\gamma },\widehat{c}^{}\left(\sqrt{F_m(\widehat{N}_c,\widehat{N}_d)}\right)^1]=1.$$
(32)
The interpretation of the relation (29) is given in Ref.?.

Our present system is also considered in terms of the boson type classical canonical variables. The canonicity condition in the TDHF theory enables us to perform this task. We introduce two sets of the variables $`(c,c^{})`$ and $`(d,d^{})`$. In the same manner as that in (I), we can show that $`c`$ and $`d`$ are related with

$$c=\gamma \sqrt{\frac{\mathrm{\Gamma }}{|\gamma |^2}\mathrm{\Gamma }^1},d=\delta \sqrt{\frac{\mathrm{\Gamma }}{|\delta |^2}\mathrm{\Gamma }^1}.$$
(33)
The expectation values of $`\widehat{N}_c`$ and $`\widehat{N}_d`$ for $`|c`$, which we denote $`N_c`$ and $`N_d`$, respectively, are given as

$`N_c`$ $`=`$ $`c|\widehat{\gamma }^{}\left(\sqrt{F_m(\widehat{N}_c,\widehat{N}_d)}\right)^1(\sqrt{m^m})^1(\sqrt{m^m})^1\left(\sqrt{F_m(\widehat{N}_c,\widehat{N}_d)}\right)^1\widehat{\gamma }|c`$ (34a)
$`=`$ $`|\gamma |^2c|(m^m)^1F_m(\widehat{N}_c,\widehat{N}_d)^1|c`$
$`=`$ $`|\gamma |^2{\displaystyle \frac{\mathrm{\Gamma }}{|\gamma |^2}}\mathrm{\Gamma }^1,`$
$`N_d`$ $`=`$ $`c|\widehat{\delta }^{}\sqrt{(m\widehat{N}_c+\widehat{N}_d+1)^1(\widehat{N}_d+1)}\sqrt{(m\widehat{N}_c+\widehat{N}_d+1)^1(\widehat{N}_d+1)}\widehat{\delta }|c`$ (34b)
$`=`$ $`|\delta |^2c|\widehat{P}_{m1}(\widehat{N}_d+1)(m\widehat{N}_c+\widehat{N}_d+1)^1\widehat{P}_{m1}|c`$
$`=`$ $`|\delta |^2{\displaystyle \frac{\mathrm{\Gamma }}{|\delta |^2}}\mathrm{\Gamma }^1.`$
The last line of each equation is directly derived from $`c|\widehat{c}^{}\widehat{c}|c`$ using (22). Thus, we have

$$N_c=c^{}c,N_d=d^{}d,$$
(35)
$`{\displaystyle \frac{\mathrm{\Gamma }}{|\gamma |^2}}\mathrm{\Gamma }^1=c|(m^m)^1F_m(\widehat{N}_c,\widehat{N}_d)^1|c,`$ (36a)
$`{\displaystyle \frac{\mathrm{\Gamma }}{|\delta |^2}}\mathrm{\Gamma }^1=c|\widehat{P}_{m1}(\widehat{N}_d+1)(m\widehat{N}_c+\widehat{N}_d+1)^1\widehat{P}_{m1}|c.`$ (36b)
The relations (33) and (35) lead us to

$`\gamma =\left[\sqrt{c|(m^m)^1F_m(\widehat{N}_c,\widehat{N}_d)^1|c}\right]^1c,`$ (37a)
$`\delta =\left[\sqrt{c|\widehat{P}_{m1}(\widehat{N}_d+1)(m\widehat{N}_c+\widehat{N}_d+1)^1\widehat{P}_{m1}|c}\right]^1d.`$ (37b)
For the relations (On the Multiboson Coherent State in Deformed Boson Scheme), we make a rough approximation :

$`\widehat{P}_{m1}|c|c,`$
$`c|G(\widehat{N}_c,\widehat{N}_d)|cG(N_c,N_d).`$ (38)
Then, we have approximate forms for $`\gamma `$ and $`\delta `$ as follows :

$`\gamma =\sqrt{m^m}\sqrt{F_m(N_c,N_d)}c,`$ (39a)
$`\delta =\sqrt{(mN_c+N_d+1)(N_d+1)^1}d.`$ (39b)
With the use of the relation (On the Multiboson Coherent State in Deformed Boson Scheme), we can give the expectation values of various operators for $`|c`$ in terms of $`(c,c^{})`$ and $`(d,d^{})`$ in the approximate forms.

The above is an outline of our idea for completing the treatment of the mutiboson coherent state given in (I). Finally, we give some remarks. In this paper, the boson operator $`(\widehat{d},\widehat{d}^{})`$ was used for the completion. Instead of the boson operator, we can use the operator $`(\stackrel{˘}{d},\stackrel{˘}{d}^{})`$ which obeys the commutation relation

$$[\stackrel{˘}{d},\stackrel{˘}{d}^{}]=1\frac{m}{(m1)!}\left(\stackrel{˘}{d}^{}\right)^{m1}\left(\stackrel{˘}{d}\right)^{m1}.$$
(40)
Under the presupposition of the existence of the state $`|0)`$ obeying $`\stackrel{˘}{d}|0)=0`$, the relation (40) leads us automatically to $`(\stackrel{˘}{d}^{})^m|0)=0`$. This means that we can formulate our idea without worrying the existence of the subspace $`\widehat{B}`$. Actually, the case $`m=2`$ was treated by the present authors. In this case, the relation (40) gives us $`\stackrel{˘}{d}\stackrel{˘}{d}^{}+\stackrel{˘}{d}^{}\stackrel{˘}{d}=1`$ and $`(\stackrel{˘}{d},\stackrel{˘}{d}^{})`$ can be regarded as fermion operator. In Ref.?, we adopted the coherent state different from the present one and, then, it is impossible to compare the results directly. The detail of the relation (40) was discussed by the present authors in relation to the Dirac canonical theory for constraint systems.