Quantum tunneling and quantum walks as algorithmic resources to solve hard K-SAT instances

We present a new quantum heuristic algorithm aimed at finding satisfying assignments for hard K-SAT instances using a continuous time quantum walk that explicitly exploits the properties of quantum tunneling. Our algorithm uses a Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_A(F)$$\end{document}HA(F) which is specifically constructed to solve a K-SAT instance F. The heuristic algorithm aims at iteratively reducing the Hamming distance between an evolving state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${|{\psi _j}\rangle }$$\end{document}|ψj⟩ and a state that represents a satisfying assignment for F. Each iteration consists on the evolution of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${|{\psi _j}\rangle }$$\end{document}|ψj⟩ (where j is the iteration number) under \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^{-iH_At}$$\end{document}e-iHAt, a measurement that collapses the superposition, a check to see if the post-measurement state satisfies F and in the case it does not, an update to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_A$$\end{document}HA for the next iteration. Operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_A$$\end{document}HA describes a continuous time quantum walk over a hypercube graph with potential barriers that makes an evolving state to scatter and mostly follow the shortest tunneling paths with the smaller potentials that lead to a state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${|{s}\rangle }$$\end{document}|s⟩ that represents a satisfying assignment for F. The potential barriers in the Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_A$$\end{document}HA are constructed through a process that does not require any previous knowledge on the satisfying assignments for the instance F. Due to the topology of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_A$$\end{document}HA each iteration is expected to reduce the Hamming distance between each post measurement state and a state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${|{s}\rangle }$$\end{document}|s⟩. If the state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${|{s}\rangle }$$\end{document}|s⟩ is not measured after n iterations (the number n of logical variables in the instance F being solved), the algorithm is restarted. Periodic measurements and quantum tunneling also give the possibility of getting out of local minima. Our numerical simulations show a success rate of 0.66 on measuring \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${|{s}\rangle }$$\end{document}|s⟩ on the first run of the algorithm (i.e., without restarting after n iterations) on thousands of 3-SAT instances of 4, 6, and 10 variables with unique satisfying assignments.

Quantum tunneling and quantum walks as algorithmic resources to solve hard K-SAT instances Ernesto Campos 1,2 , Salvador E. Venegas-Andraca 1* & Marco Lanzagorta 3 We present a new quantum heuristic algorithm aimed at finding satisfying assignments for hard K-SAT instances using a continuous time quantum walk that explicitly exploits the properties of quantum tunneling. Our algorithm uses a Hamiltonian H A (F) which is specifically constructed to solve a K-SAT instance F. The heuristic algorithm aims at iteratively reducing the Hamming distance between an evolving state |ψ j � and a state that represents a satisfying assignment for F. Each iteration consists on the evolution of |ψ j � (where j is the iteration number) under e −iH A t , a measurement that collapses the superposition, a check to see if the post-measurement state satisfies F and in the case it does not, an update to H A for the next iteration. Operator H A describes a continuous time quantum walk over a hypercube graph with potential barriers that makes an evolving state to scatter and mostly follow the shortest tunneling paths with the smaller potentials that lead to a state |s� that represents a satisfying assignment for F. The potential barriers in the Hamiltonian H A are constructed through a process that does not require any previous knowledge on the satisfying assignments for the instance F. Due to the topology of H A each iteration is expected to reduce the Hamming distance between each post measurement state and a state |s� . If the state |s� is not measured after n iterations (the number n of logical variables in the instance F being solved), the algorithm is restarted. Periodic measurements and quantum tunneling also give the possibility of getting out of local minima. Our numerical simulations show a success rate of 0.66 on measuring |s� on the first run of the algorithm (i.e., without restarting after n iterations) on thousands of 3-SAT instances of 4, 6, and 10 variables with unique satisfying assignments.
K-SAT is a most important problem in computer science which is known to be NP − Complete for all K > 2 1,2 . We focus on the set of satisfiable instances of the K-SAT, i.e. those instances for which there is at least one assignment that returns a logical 1. Satisfying a K-SAT instance may be easy or require very significant efforts, it all depends on its number of variables and clauses. Instances with one or few satisfactory assignments are harder to solve because satisfactory assignments live in a set with cardinality 2 n , where n is the number of binary variables upon which the instance is built. Those K-SAT instances for which there is only one satisfying assignment constitute the focus of this paper.
A central activity in quantum computing consists of using quantum mechanical properties as resources in order to build algorithms that may outperform their classical counterparts. Along these lines, we find several quantum algorithms that have been designed to solve K-SAT instances [3][4][5] . In this paper, we present a quantum algorithm in which continuous quantum walks and quantum tunneling are explicitly used as computational resources.
Originally designed to model quantum phenomena [6][7][8][9] , quantum walks are an advanced tool for building quantum algorithms (e.g. [10][11][12][13][14][15][16][17][18] ) that has been shown to constitute a universal model of quantum computation 19,20 . Quantum walks come in two variants: discrete and continuous 21 . The continuous-time quantum walk was proposed by Childs et al 11 and the continuous time quantum walk on a hypercube was first studied by Kendon and Tregenna 22 and later extensively analyzed by Drezgich et al 23 .
Quantum tunneling is a quantum mechanical phenomenon in which a particle passes through a potential barrier. The probability of finding a particle inside the potential barrier decreases exponentially, which in turn increases the probability of finding the particle outside the potential 24 , as if the particle were expelled from the barrier. This behavior suggests that quantum tunneling can be used as an algorithmic tool in which researchers

K-SAT problem
The K-SAT problem is one of the most important problems in computer science. It is derived from the SAT problem, the first problem proved to be NP − complete 1 . K-SAT is a constraint satisfaction problem where a Boolean formula has to be satisfied, i.e. its result has to be true. A K-SAT instance is written in conjunctive normal form (CNF), that is, a conjunction of clauses, where a clause is a disjunction of literals, and a literal can be a logical variable or its negation.
Formally speaking, the K-SAT problem can be stated as follows: given a K-SAT instance P, is there an satisfying truth assignment for P? In other words, the K-SAT problem consists of finding assignments that satisfy a given CNF Boolean formula, where K is the number of literals in each clause. For example: where ∨ is the logical symbol for OR, ∧ is the logical symbol for AND, {x 1 , . . . x 6 } is the set of logical binary variables, and x j is the negation of x j . This is a 3-SAT instance since each clause contains exactly 3 variables, and the binary string x s = [1, 1, 0, 0, 1, 0] is one of the assignments that satisfy P(x) . This can be verified by substituting x s into P(x) Satisfying a K-SAT instance P can be easy, difficult or impossible, depending on the number of bit strings that make P true. For example, let us suppose we have a K-SAT instance P defined over n binary variables. If we run a brute-force approach on P, we may simply substitute each and every bit string from x 0 = 000 . . . 0 to x 2 n −1 = 111 . . . 1 on P, only stopping either when finding P(x s ) = 1 or when we have substituted all bit strings from x 0 to x 2 n −1 without finding any satisfying assignment. If P has many satisfying assignments, we will eventually find one of those assignments, possibly sooner than later. However, if P has only a few satisfying assignments (along these lines, the worst case scenario would be having only one satisfying assignment), finding them will eventually happen but it may well take a substantial (i.e. exponential) amount of time. If no bit string from x 0 to x 2 n −1 satisfies P, then P is unsatisfiable but it took us an exponential amount of resources to find that out.
As stated above, the K-SAT is an element of the set of NP − Complete problems, which in turn is a subset of NP (the abbreviation of "Non-deterministic Polynomial Time"), the set of problems in which a solution can be tested in polynomial time. Being an NP − Complete problem means that an arbitrary instance of any problem under the NP category can be rewritten as an instance of an NP − Complete problem in a polynomial number of steps 2,33,34 . The first problem proved to be NP − Complete was the SAT problem, and later the K-SAT, which is a variation of SAT, was proven to be also NP − Complete for K > 2 1,2,34-36 . If an algorithm that runs on a Deterministic Turing Machine in polynomial time is discovered for an NP − Complete problem, it implies every NP problem can be solved in polynomial time, hence proving P = NP.

Continuous time quantum walk on a hypercube graph
An n-dimensional hypercube graph can be constructed by assigning each vertex to one of the 2 n possible n-variable binary combinations, and connecting the nodes that have a distance of 1 in the Hamming distance (the number of positions/bits at which the corresponding elements between two arrays are different). The Hamiltonian H A used in our algorithm describes the dynamics of a quantum walk over a hypercube with potential barriers. The reason why we chose to use the hypercube graph is to reduce the distance between |ψ j � and |s� with every iteration of the algorithm by taking advantage of the hypercube topology, where closer states have a shorter Hamming distance. During the evolution of an n-qubit state |ψ j � (j is the iteration number) under e −iH A t , most of the probability flux will tunnel through the paths from |ψ j � to |s� , and after a time t f the system is measured and the superposition collapses. If the post-measurement state |ψ j+1 � is different from |s� , there is a good chance |ψ j+1 � is a state in a tunneling path from |ψ j � and |s� which thanks to the hypercube topology will be at shorter Hamming distance to |s� , h H (|ψ j+1 �, |s�) < h H (|ψ j �, |s�).
The structure of the the adjacency matrix A(n) for an n-dimensional hypercube can be written as If a state that is the tensor product of n qubits in the computational basis evolves under the unitary operator e −iA(n)t (which is equivalent to evolving under H A without potential barriers), the probability at some time of measuring a state |z� , that is also an n-qubit tensor product, can be calculated as follows. Firstly, we need to realize that the Hermitian matrices A i commute with each other, for example which means we can express e −iA(n)t as a product of unitary operators e −iA(n)t = e −iA 1 t e −iA 2 t ...e −iA n t . Then, by expanding e −iA j t into a Taylor series and gathering the terms, we obtain which means e −iA j t is equivalent to only evolving the ith qubit under e −iσ x t . As a result, we get e −iA(n)t = e −iA 1 t e −iA 2 t ...e −iA n t = n j=1 e −iσ x t . The probability P z (t) of measuring a state |z� from |x� evolving under e −iA(n)t depends only on the Hamming distance between |x� and |z� . If we evolve |x� = |0� ⊗n , the probability of finding an arbitrary state |z� is given by where n 0 and n 1 are the number of 0's and 1s appearing in |z� respectively, meaning n 1 is also the Hamming distance between |x� and |z�.

Formal description of our algorithm
Our algorithm is formally introduced in ALGORITHM 1. A key element of our algorithm is the construction of the Hamiltonian H A (F, |ψ j �, j) , which determines the dynamics under which a state |ψ j � will evolve during the jth iteration of the algorithm. Each iteration of the algorithm consists of the evolution of a state |ψ j � during a time t f (this is the measurement frequency, for the cases analyzed in this paper t f = 3 2 π as discussed in "Numerical calculation of parameters and additional simulations") followed by a measurement that collapses the superposition, and then a check of the post measurement state |ψ j+1 � to see if it represents a satisfying assignment for F (the checking can be done in O(mK)), in case it does not, H A gets updated for the next iteration.
The Hamiltonian H A for the jth iteration is where A(n) is the adjacency matrix of an n-dimensional hypercube graph and n is the number of variables in the instance F, which can be constructed using Eq. (1), and is the potential part of H A . The components of V are: • |ψ j � is the tensor product of n qubits in the computational basis. |ψ j � represents an assignment ψ j for the instance F. |ψ 0 � is an initial random state, while |ψ j � for j > 0 is the post measurement state for the jth iteration. • α = α(F, ψ j ) is the number of clauses in F that are not satisfied by the assignment represented by |ψ j �.
• M = M(F) is a diagonal matrix that has as the element corresponding to |a��a| the number of clauses in F not satisfied by the assignment represented by |a� . It is worth noting that the elements of M(F) arise naturally from a construction process (shown later in this section) without any prior knowledge on the satisfying solutions of F. • γ (n, m, K) and δ(n, m, K) are both positive real parameters dependent in the number of variables n, the number of clauses m, and the number of literals per clause K in the instance F. The reason of their dependency on n, m, and K is shown later in this section. γ increases an amount δ every iteration which, as we explain later in this section, helps with the amplification of the state |s� that represents a satisfying assignment for F for an evolution during a time t f (the measurement frequency). We do not have an analytical expression for γ and www.nature.com/scientificreports/ δ , so we have numerically estimated their values for some cases. The values and their effects are discussed in more detail in "Effective Hamiltonian via degenerate perturbation theory" and "Numerical calculation of parameters and additional simulations".
The rationale behind the elements of H A is the following. The potential V of Hamiltonian H A is a diagonal matrix with potential values specific to each state. The potential barriers between the states are designed to increase the tunneling between states |ψ j � and |s� . The reason most of the tunneling occurs between |ψ j � and |s� is because the elements in the potential V corresponding to the outer products |ψ j ��ψ j | and |s��s| are the only ones with a value of 0, which is a consequence of the construction process of V.
To construct the potential V, the first step is to build the Non-Satisfied Clauses matrix M(F). As stated before, M(F) is a diagonal matrix whose non-zero entries are the number of clauses in F which an assignment a does not satisfy. Assignment a is represented by state |a� and its corresponding non-zero entry in M(F) is located on the position that corresponds to the outer product |a��a|.
To s h o w h o w M ( F ) i s b u i l t , w e n o w i n t r o d u c e a n e x a m p l e b a s e d o n , and the assignments that do not satisfy C 1 are : The states that represent the assignments that do not satisfy C 1 are |000� and |001� , and the sum of their outer products is The NS(C 1 ) matrix form is given by diag(1, 1, 0, 0, 0, 0, 0, 0) . Note that we can use I when doing the outer product summation to represent the variables that do not appear in a clause ( x 3 in the case of C 1 ), hence reducing the number of elements needed to construct each NS(C i ) . The number of non I elements in every NS(C i ) is only dependent on K, the number of binary variables on each and every clause of a K-SAT instance, which is fixed for a K-SAT instance. This allows us to construct each NS(C i ) with a non-exponential number of elements (see "Analytical approximation of the algorithm behavior" for a detailed explanation of how to efficiently construct NS(C i ) matrices). NS(C 1 ) is a diagonal matrix with 1s in the elements corresponding to the outer products of the states that do not satisfy C 1 . In this case the 1s are clearly in the positions corresponding to the outer products of |000� and |001� , while the zeroes correspond to the outer products of the states that represent a satisfying assignment for C 1 .
If the matrices NS(C i ) are calculated for every clause in an instance F, the diagonal elements that are 0 in every NS(C i ) correspond to the outer products of the states that satisfy every clause, thus satisfying instance F. The sum of every NS(C i ) results in the diagonal matrix where an entry in M(F) corresponding to an outer product |a��a| is equal to the number of clauses not satisfied in F by the assignment represented by |a� , and m is the number of clauses in F. For our example M(F 1 ) is which in matrix notation is given by From Eq. (7) we can see that the assignment corresponding to |000� does not satisfy 1 clause in F 1 , |001� does not satisfy 1 clause, |010� does not satisfy 2 clauses, and so on. There is only one zero in the diagonal of M(F 1 ) , which corresponds to the outer product of |111� , meaning the assignment Since the probability of tunneling is higher between states at potentials with similar values, we want to start with an initial state in a very small potential, expecting its probability flux to mostly tunnel through the shortest tunneling paths with the smaller potentials leading to |s� that represent a satisfying assignment for an instance F (the element corresponding to |s��s| in M(F) is 0, the global minimum) thus amplifying the probability of measuring |s� . Since we do not know a priori which states correspond to a small potential, we can just set to 0 the potential corresponding to the outer product |ψ j ��ψ j | by classically calculating the number of clauses α(ψ j , F) the assignment corresponding to |ψ j � does not satisfy (this can be done in O(mK)), and then subtracting the outer product |ψ j ��ψ j | multiplied by α(ψ j , F) from the matrix M(F). As a result we get a 0 in the element corresponding to |ψ j ��ψ j | in M(F).

From Eq. (3), the Hamiltonian H A is
Suppose |ψ j � = |000� . We build A(n) using Eq. (1). Then, the resulting H A is given by www.nature.com/scientificreports/ Note that explicitly calculating every element in M(F) consumes an exponential amount of resources in the number of qubits. We do it explicitly in this paper with the purpose of showing the logic and behavior of the algorithm. In a real world implementation, M(F) would be a physical process that acts as a blackbox.
Parameters discussion. As stated before, γ and δ are both real positive parameters. First we will explain the effect γ + δj has in the probability of measuring |s� , then we proceed to explain why it increases with every iteration, and finally why they are dependent on n, m and K, the number of variables, clauses, and literals per clause in F respectively.
• The effect γ + δj has on H A is clearly a change in the size of the potential barriers. If (γ + δj) < 1 the size of the potential barriers get reduced which makes a state |ψ j � evolving under e −iH A t to tunnel and amplify |s� faster than if (γ + δj) = 1 . This occurs because the smaller the potential barriers the easier it is for the probability flux to transmit through them, but it has the unwanted side effect of letting more of the probability flux to transmit to states outside the shortest paths from |ψ j � to |s� which translates into a smaller amplification of |s� . On the other hand, if (γ + δj) > 1 the time it will take |ψ j � to tunnel to |s� increases but has the benefit of less probability flux getting out of the main tunneling paths from |ψ j � to |s� which in turn increases the amplification and also reduces the small oscillations created by the reflections in the potentials and the probability flux transmitting out of the main tunneling paths. Simulations and analytical approximations showing the previously mentioned effects are shown in "Numerical calculation of parameters and additional simulations". • To understand why parameter γ gets an increment δ on every iteration, it is required to understand how the time it takes to get the maximum amplification of |s� varies depending on its Hamming distance to |ψ j � , d H (|s�, |ψ j �) . Intuitively, a state |ψ j � evolving under e −iH A t will amplify faster |s� the smaller d H (|ψ j �, |s�) is, since it has to tunnel through less potential barriers (see "Numerical calculation of parameters and additional simulations" for simulations and approximations exhibiting this behavior). As explained above, thanks to the hypercube topology we have a good probability of getting a post measurement state |ψ j+1 � from the tunneling paths from |ψ j � to |s� with a smaller Hamming distance to |s� compared to |ψ j � , d H (|ψ j+1 �, |s�) < d H (|ψ j �, |s�) . So, the purpose of increments δ to γ on every iteration is to keep approximately the same value for the measurement frequency t f , so that the probability P s (t) of measuring |s� is greatly amplified for multiple distances d H (|ψ j �, |s�) approximately at the same time. • Parameters γ , δ depend on n, m and K (the number of variables, clauses and literals per clause of F respectively) because as n increases, the average number of iterations also increases. This is because the average Hamming distance is also increased, and having greater γ and δ can make tunneling too difficult from longer Hamming distances in few iterations. As for the dependency on m and K it is because we cannot know the value of every region in the potential V since they are 2 n values, but the average potential value apv = m 2 K is only dependent on m and K for (γ + δj) = 1 (the derivation of apv is presented in "Numerical calculation of parameters and additional simulations", Eq. (11)). If an instance F 1 has a greater apv than an instance F 2 , and both instances have the same number of variables, we can expect the optimal values of γ and δ for F 1 to be smaller than the optimal values for F 2 since smaller values of γ and δ will facilitate the tunneling through the bigger potential barriers found in H A (F 1 ) . We have estimated, via digital computer algorithms, optimal values of γ and δ , by doing a sweep over the values, for 3-SAT instances with 4, 6 and 10 variables with 16, 24, and 40 clauses respectively. Table 1 shows smaller values for larger instances with more variables in accordance to results presented above. The details on the numerical approximations of γ and δ are shown in "Conclusions".
We have also numerically estimated, for different values of γ and δ , the time during which we evolve the system t f , the measurement frequency, giving the same value t f = 4.72 in every case (details in "Conclusions"). www.nature.com/scientificreports/ The analytical approximation of t f via perturbation theory (shown in "Numerical calculation of parameters and additional simulations") gave t f = 3 2 π , in close accordance with the numerical results.

Implementation details
The algorithm may be efficiently implemented by approximating the evolution operator as a set of quantum gates by using techniques like the Trotter-Suzuki decomposition 37 or its variations. We have developed this technique to show that there is at least one efficient method to implement our algorithm. It remains, as future work, to identify optimal implementation methods that may preserve the continuous nature of the algorithm 38 .
H A can be decomposed as a sum of Hermitian matrices acting on small subspaces where NS(C i ) are the Hermitian matrices that compose M(F). M(F) and A(n) (Eq. (1)) seem to be expensive to construct, but in reality M(F) can be constructed using O(m) unitary matrices, while A(n) can be constructed in O(n). Each NS(C i ) can be built from 2 K unitary matrices acting on a small subspace, this is because K is the number of variables contained in each and every clause of a K-SAT instance and, consequently, 2 K is a rather small number, certainly much less than 2 n , where n is the number of binary variables over which a K-SAT instance is defined. So, for any instance of the K-SAT problem K is a fixed integer number equal to the number of variables in each and every clause, then 2 K = O(1) for any given K-SAT instance. For example, the first clause of ) and the states that represent the assignments that do not satisfy C 1 are |000� and |001� . The sum of their outer products is the matrix NS(C 1 ) where σ z is the Pauli operator σ z = |0��0| − |1��1| . So, M(F) can be expressed as the sum of 2 K m Hermitian operators. In general, the number of Is in NS(C i ) is n − K , in our example (Eq. (9)) K = 2 and n = 3 so it has one I. For instances with more variables, the number of I increases which means that each NS(C i ) only acts in a small subspace. Expressing NS(C i ) as the sum of 2 K Hermitian operators, as in Eq. (10), requires using K2 K−1 Pauli operators σ z . The number of σ z operators in NS(C i ) is independent from the number of variables in the instance, if the instance had more variables it would just add more Is to each NS(C i ) . Consequently, we do not need an exponential number of Pauli matrices to construct M(F).
The unitary e iα|ψ j ��ψ j |t can be efficiently implemented using ancilla qubits. A controlled gate applies a NOT gate to an ancilla qubit initially in |0� if the original n qubits are in |ψ j � (this can be done with a polynomial number of gates and ancillas 39 ). Then, if the ancilla qubit is in |1� a controlled gate applies e iαt to any of the original n qubits. The required projector for lowering the energy of |ψ j � can also be built with a gadget proposed by Dodds et al 40 .

Analytical approximation of the algorithm behavior
To analytically approximate the probability P s (t) of finding a satisfying assignment to a K-SAT instance F, which is equivalent to having |s� as post-measurement quantum state, and to understand how γ + δj affects this probability, we use the effective Hamiltonian H eff obtained by degenerate perturbation theory. The effective Hamiltonian only acts on a reduced space and only describes part of the energy spectrum of the complete Hamiltonian. It has the advantage of being simpler than the complete Hamiltonian and it is more general than normal perturbation theory since it takes into account the effects of the off-diagonal elements that link the states of interest. For example, it is used in the Born-Oppenheimer approximation to separate the electronic wavefunction from the nuclear wavefunction of a molecule 41 . In our case we want to separate the subspace spanned by |ψ j � and |s� in H A to later approximate P s (t).
Perturbation theory is a set of approximations used to describe a complex quantum system in terms of a simpler one. The approximation is based on the idea of using a Hamiltonian (the unperturbed Hamiltonian H 0 ) with known solutions and adding a perturbation Hamiltonian H 1 that slightly modifies the system. The solutions of the perturbed Hamiltonian H p can be expressed as corrections to the solutions of the unperturbed Hamiltonian H 0 . Degenerate perturbation theory is a variation of perturbation theory used for unperturbed Hamiltonians with degeneracy in their energy spectrum.
The perturbed Hamiltonian is H p = H 0 + H 1 where is a dimensionless parameter small enough so that the spectrum of H 0 constitutes a good starting point to approximate the spectrum of H p . In the following examples is set to 1 (as it is frequently done in the literature) in the understanding that the elements of H 1 are smaller than those of H 0 .
As for our algorithm, to approximate the Hamiltonian already presented in Eq. (3) We have chosen H 0 and H 1 to be V and A(n) respectively because the entries of H 1 = A are much smaller than those of H 0 = V . Entries of H 1 = A are composed by 1s and 0s while entries of H 0 = V have average value potential value calculated as follows: For big instances, the contribution of −α|ψ j ��ψ j | can be ignored since it only affects one of the 2 n possible states. So, Also, from Eq. (6), we know that So, the average potential value is the trace of M(F) divided by 2 n .
To compute tr(M(F)) we need to know how many elements in M(F) each of the m NS(C i ) matrices contribute ( C i are the clauses in F). As it can be seen in "Implementation details" Eq. (5), the trace of each NS(C i ) is tr(NS(C i )) = 2 n−K , where K is the number of literals per clause, so which means that the larger the number of variables and clauses, and the higher value of γ + δj , the better the approximation.
It is convenient to use H 0 = V since it is a diagonal matrix and calculating its energies and eigenstates is straightforward. Also, there is degeneracy in the energy spectrum of V since the energies of |ψ j � and |s� (both eigenstates of H 0 = V ) are E s = E ψ j = 0.
Effective Hamiltonian via degenerate perturbation theory.. The effective Hamiltonian H eff acts on a reduced space and only describes part of the energy spectrum of the complete Hamiltonian, which in our case is H A . The effective Hamiltonian we want to approximate is the one spanned by the states |ψ j � and |s� since H A has two eigenvectors that can be approximated as a linear combination of |ψ j � and |s� by the use of the effective Hamiltonian H eff . The approximated eigenvectors and energies then will be used to approximate P s (t) , the probability of getting |s� as post-measurement state, i.e., finding a satisfying assignment for a K-SAT instance F. As before, our analysis is focused on instances with a unique satisfying solution to make it clearer and because those instances are the hardest to solve.
States |ψ j � and |s� may be far from each other in terms of the Hamming distance. If the Hamming distance between |ψ j � and |s� is d, i.e. d H (|ψ j �, |s� = d , then a dth order approximation to H eff is needed for the effects of quantum tunneling between |ψ j � and |s� to be present in H eff . In this analysis we use a 3rd order approximation but the process for higher orders can be easily generalized.
From Gottfried 42 , the spectrum of H 0 = V has degeneracy from the energies E s = E ψ j = 0 corresponding to the states |ψ j � and |s� . We will call D the degenerate subspace in H 0 = V spanned by the states {|ψ j �, |s�} ∈ D.
The perturbed states |a� that correspond to the degenerate states |α� (where |α� represent any of the states in D) that split due to the perturbation can be written as the following linear combination where |µ� are the states outside D.
To obtain the effective Hamiltonian we will approximate d µ and do some manipulations so we end up with an eigenvalue problem in the subspace D. As we will see, the order of coefficients c α are O(1) while coefficients d µ are O( ) . Based on this, we will approximate the states |a� as a linear combination of the degenerate states |α�.
Using   (14) is O( 2 ) and can be drooped for the approximation. Since the difference between E α and E a is negligible compared to E ν − E a , E a can be replaced by E α = 0 Substituting (15) in (14) we obtain a better approximation of d ν than (15). We can iteratively repeat this process to generate higher order approximations of d ν (and in turn, as we will see in Eq. (18), higher order approximations of H eff ).
Noting d ν and d µ both represent the coefficients of the states outside D, we can substitute (16) in (13) This is equivalent to the eigenvalue problem (13) For the third order approximation of the effective Hamiltonian H Approximation of the probability of finding a satisfying assignment to a K-SAT instance. We now present an approximation of the probability P s (t) , the probability of finding a satisfying assignment to a K-SAT instance F, which is equivalent to obtain |s� as post-measurement state of a quantum evolution process starting with the quantum state |ψ j � . For this purpose, we use the effective Hamiltonian H eff from Eq. (18) for three different cases: d H (|ψ j �, |s�) = 3 , d H (|ψ j �, |s�) = 2 , and d H (|ψ j �, |s�) = 1 . The approximation for P s (t) for distances not presented here can be easily obtained by slightly modifying the procedure.
First we show the approximation H where |µ� and |ν� are eigenvectors with eigenvalues E µ , E ν in H 0 = V different from |ψ j � and |s� . The resulting effective Hamiltonian will be 2 × 2 since it is from the subspace spanned by {|ψ j �, |s�} from which we will approximate two eigenvectors of H A as a linear superposition of |ψ j � and |s� , and later use them to approximate P S (t).
The linear term in (19) has no contribution to the diagonal of H eff since A (as any adjacency matrix) has a diagonal full of 0s. As for the off diagonal elements of the linear term, there is no way for |ψ j � and |s� to tunnel from one to the other without going through intermediate states in the case d H = 3 , so there are no off diagonal elements contributed by the lineal term.
Evaluating the second order term, we can see it only has diagonal elements, the contributions to the diagonal come from evaluating the term with states at d H (|µ�, |s�) = 1 or d H (|µ�, |ψ j �) = 1 . Each contribution to this term is divided by its corresponding E µ which is a product of γ + δj and the number of clauses the corresponding state |µ� does not satisfy. www.nature.com/scientificreports/ The third order term is the only contributor to the off-diagonal elements of H eff in the case d H = 3 . These contributions come from evaluating the third order term with the states that form the shortest tunneling paths from |ψ j � to |s� , these states are the ones at a distance d H = 2 from |ψ j � or |s� and d H = 1 from the other. The easier it is for the states |ψ j � and |s� to tunnel through these paths (when the potential barriers are smaller) the greater the value of the off-diagonal elements will be. We can check the previous statement from Eq. (19) since every contribution from each tunneling path is 1 divided by E µ E ν .
The where k ± and k ′ ± the are normalization constants and can be related as follows: The approximated eigenstates of H A , |φ ± � written as a linear combination of |ψ j � and |s� are From Eq. (24) and orthogonality, the eigenvectors can be rewritten as where g is a normalization constant. We can use |φ ± � to approximate the evolution of the real eigenstates of the Hamiltonian H A . The evolution of a true eigenstate of H A is by substituting the eigenstates and energies by their approximation we obtain Using it, we can get an approximation of probability P

That is,
The frequency of the probability P It is easy to follow a similar process to calculate P (2) s (t) , the case with d H (|ψ j �, |s�) = 2 . In this case at least a second order approximation is needed for the effects of tunneling (the off diagonal elements) to appear in H (2) eff .
The corresponding H eff is slightly different to the previous case which has an impact in the frequency of P s (t) . The first order contribution to H eff , as in the previous case, is also null for the same reasons. In this case the second order term gives all the elements for H eff . The contributions from the second order come from evaluating the states at d H (|µ�, |ψ j �) = 1 or d H (|µ�, |s�) = 1 . Some of these states will be part of the shortest tunneling paths from |ψ 0 � to |s� , meaning d H (|µ�, |ψ j �) = d H (|µ�, |s�) = 1 , and are the ones that contribute to the off-diagonal elements. As in the previous case ( d H = 3) , the easier it is for |ψ j � to tunnel through a barrier to get to |s� the greater the contribution from that path. Each contribution is divided by the corresponding E µ , that is a product of γ + δj and the number of clauses the corresponding state |µ� does not satisfy. This gives the following effective Hamiltonian Following a process similar to the one presented for d H = 3 , we obtain s (t) frequency, when γ + δj >> 1 , remembering � = 1 γ +δj , P s (t) has the term 4v 2 12 inside the square root compared to P s (t) where 2 4v 2 12 vanishes. Also, the off-diagonal elements v 12 from H eff when d H = 2 tend to be greater compared to the case d H = 3 since there are less potential barriers in the tunneling path from |ψ j � to |s� . Consequently, in most cases the frequency of P (2) s (t) will be higher than he frequency of P   (27) does not have dependency on and has a period of π . Hence, for t = π 2 + πm , where m ∈ Z , P (1) s (t) is at its highest value. As γ + δj increases, the potential barriers around the states |ψ j � and |s� get larger, reducing the tunneling to the states outside of D (for d H = 1 the states in D are neighbors) which makes the real behavior of P   With this in mind, the values of γ and δ have to keep the evolution time t f = π 2 + πm , where m ∈ Z , as a time that greatly amplifies P s (t) for multiple values of d H (|ψ j �, |s�) . In fact, our numerical simulations using the optimal parameters, found by doing a sweep, show this property of preserving t f as an evolution time that greatly amplifies P s (t) . Our simulations give us an optimal measurement frequency of t f = 4.72 , that is approximately 3 2 π , in agreement with analytical results. Keep in mind that for larger numbers of variables and clauses this value may change to t f = π 2 + πm , with a different value of m.
Comparisons between the analytical approximation and simulations. In this section, we compare the analytical approximations of P s (t) from the previous section to our numerical simulations of P s (t) for the distances d H (|ψ j �, |s�) = 1, 2, 3 . The instance used for these comparisons is the following 6 variable instance with 24 clauses and the unique satisfying assignment s = [110111].
To exemplify the case d H = 3 we used the state |ψ j � = |010001� . The effective Hamiltonian using degenerate perturbation theory up to the third order, as in Eq. (20), is Figure 1 compares analytical approximation and numerical simulations of P s (t) evolving from |ψ j � = |010001� for various values of γ + δj . As the value of γ + δj increases the approximation gets better since the tunneling www.nature.com/scientificreports/ to states outside the main tunneling paths from |ψ j � to |s� decreases, and the reflections at the potential barriers also decrease. This contributes to reduce the small oscillations seen in the numerical simulations of P s (t) . We can also see that even for small values of γ + δj , the approximation of the period is accurate. For d H = 2 , we use the initial state |ψ j � = |110001� . By using perturbation theory up to the second order, as in Eq. (26), we obtain the following effective Hamiltonian Figure 2 shows comparisons between the analytical approximation and the numerical simulations of P s (t) evolving from |110001� for various values of γ + δj.
Lastly, for d H = 1 and using the initial state |110011� the effective matrix is simply Figure 3 shows comparisons between the analytical approximation and the numerical simulations of P s (t) evolving from |110011� for various values of γ + δj Graphs from Figs. (1,2,3) show the advantage of including potential barriers compared to the exponentially small overlap P s (t) = sin 6 (t) cos 6 (t) = sin 6 (2t)/2 6 , obtained from Eq. (2), corresponding to solely using the hypercube adjacency matrix as the quantum walk Hamiltonian.  In general, the period of P s (t) grows longer as d H increases as a consequence of having to tunnel through longer paths, which is also reflected in the analytical approximation for the period by becoming more dependent on γ + δj as a consequence of H eff having the mathematical structure where d is the distance d H (|ψ j �, |s�) = d and also the minimum approximation order necessary for the tunneling effects to be present in the approximation. As d H = d grows the off-diagonal elements in H eff that serve as couplings between |ψ j � and |s� get smaller which implies longer periods for transitioning between states.

Numerical calculation of parameters and additional simulations
So far, we have been able to numerically estimate optimal values for parameters γ and δ ; moreover, our analytical approximation of t f = π 2 + πm , m ∈ Z has multiple possible values. Since it remains as future work to derive analytical formulae for these parameters, we have approximated three sets of these parameters for 3-SAT instances of 4, 6, and 10 variables, with 16, 24, and 40 clauses respectively.
Our database, fully available in 32 , consists of thousands of 3-SAT instances from which 5799 have 4 variables and 16 clauses, 1949 have 6 variables and 24 clauses, and 606 have 10 variables and 40 clauses. Every instance Figure 2. Comparison between simulations (blue line) and the analytical approximation (orange) of P s (t) for initial state |110001� , d H = 2. www.nature.com/scientificreports/ in the database has a unique satisfying solution. All instances were randomly created by a uniform distribution, and then solved by exhaustive substitution as our instances are small enough for that. Instances with unique satisfying assignment were saved for the database. The sets of parameters were computed by doing thousands of simulations on 3-SAT instances for the three different number of variables in our data set. The parameters were computed by slightly varying the value of each parameter and running 500 simulations for each combination of parameters. Half the dataset was used to find the parameters while the other half was used to test their performance. Each simulation was run using a randomly selected instance from the dataset, starting with a randomly selected initial state, and stopped after n measurements, where n is the number of variables of the instance, meaning only the first run of the algorithm (without resetting after n iterations) was considered for these calculations.
In every case the measurement frequency t f was around 4.72 which is approximately 3 2 π in accordance to our analytical approximation. Future versions of the algorithm may make use of variable measurement times, possibly by the use of techniques used in the literature of first detection time for quantum walks 43,44 .
The optimal sets of parameters we estimated for each number of variables and clauses are shown in Table 2. Alongside them we have included statistical information from 1000 simulations for each of the three different number of variables: The success rate of having |s� as post-measurement quantum state (that is, the success rate of finding the unique solution of each instance used in this study) for the first run of the algorithm (n iterations at most), average number of necessary iterations to get |s� as post-measurement quantum state, and the standard deviation in the number of iterations to get |s� as post-measurement quantum state.
Taking in account all the simulations we ran for Table 2, the success rate was 0.66 for the first run of the algorithm. The success rate for the 10 variable instances may seem low, but remember it is only for the first run, www.nature.com/scientificreports/ for the second run its success rate increases to 0.63 and the fifth run is already 0.93. Also, we are using linear increments for γ + δj , the use of a different function may translate into a better performance. Figure 4 presents simulations of P s (t) for two randomly selected 10 variables 3-SAT instances from the database 32 showing the expected behavior for runs where the satisfactory state |s� is measured.

Conclusions
We have presented a quantum walk-based algorithm that explicitly utilizes potential barriers and quantum tunneling to find satisfying assignments for hard K-SAT instances.
One of the novelties of our paper is to harness by design quantum tunneling as a computational resource and to use it as a controllable feature that can be engineered to provide additional processing power to a quantum algorithm. Furthermore, we have derived a detailed mathematical procedure to build a family of Hamiltonians to solve K-SAT instances via continuous quantum walks with potential barriers and quantum tunneling. Based on the structure of these Hamiltonians, we have used degenerate perturbation theory to derive an analytical approximation of the probability of finding satisfying assignments to K-SAT instances under our algorithmic approach.
As for the set of parameters required to run our algorithm, we have made an extensive numerical study that has allowed us to present good estimates of those parameters. It remains as future work to derive analytical formulae for this set of parameters, a future work which is encouraged by our numerical results.  www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.