The QAOA aims to find good approximate solutions to combinatorial optimization problems (Farhi et al., 2014) that are defined by a classical objective function E(x) mapping N-binary strings to real values: E x : 0,1 N → R . (1) The algorithm’s goal is to find a binary string x* that minimizes the function E _min ≡ min _x E(x), or at least that achieves a good approximation ratio r ∈ [0, 1]. max x E x − E x max x E x − min x E x ≥ r . (2) This optimization problem is equivalent to finding the ground state of a spin Hamiltonian E ̂ , where each binary variable x _i describes the state of one spin σ i z as x → 1 2 ( 1 + σ z ) : E ̂ min ≡ min σ z 〈 σ z | E ̂ σ z | σ z 〉 = min x 〈 x | E ̂ x | x 〉 . (3)

A well-known approach to solving these problems is the quantum adiabatic algorithm (QAA), which ensures the achievement of the global minimum given a sufficiently long run time T (Farhi et al., 2000; Farhi et al., 2001). The adiabatic protocol is guided by a time-dependent Hamiltonian, such as H ̂ QAA t = − 1 − f t T H ̂ x + f t T E ̂ , (4) where H ̂ x = ∑ i N σ i x with σ i x being the Pauli operator acting on the ith qubit and f(x) is the schedule function with f(0) = 0 and f(1) = 1. The adiabatic theorem guarantees that a sufficiently slow time evolution will keep the system in its ground state. Thus, an evolution that starts with the ground state of H ̂ x , | Ψ 0 〉 = 1 2 N ∑ x | x 〉 ≡ | + 〉 ⊗ N , (5) will approximately bring the system to the ground state of E ̂ , | Ψ T 〉 = | x min 〉 . (6) To ensure success, the runtime T must scale as T = O Δ min − 2 , where Δ_min is the minimum spectral gap (Albash and Lidar, 2018). The QAOA proposes a Trotterized approximation to adiabatic evolution consisting of a quantum circuit built by the alternation of the following two operators: | Ψ QAOA - p γ , θ 〉 = U H ̂ x , θ p U E ̂ , γ p … U H ̂ x , θ 2 U E ̂ , γ 2 U H ̂ x , θ 1 U E ̂ , γ 1 H ⊗ N | 0 〉 ⊗ N , (7) with H denoting single-qubit Hadamard gates, and U E ̂ , γ = e − i γ E ̂ , (8) U H ̂ x , θ = e − i θ H ̂ x = ∏ j N e − i θ σ j x ≡ R x θ ⊗ N , (9) where γ = (γ ₁, γ ₂, …, γ _p ) and θ = (θ ₁, θ ₂, …, θ _p ) are set of variational angles that can be tuned to approximate the ground state of the cost Hamiltonian E ̂ and p is the number of layers that determines the depth of the quantum circuit and the accuracy of the algorithm. The infinite depth limit p → ∞ returns the adiabatic evolution, showing that we can get a good enough approximation to the minimum of the optimization problem for sufficiently large p.

Let us highlight the QAOA’s potential in small-depth regimes by introducing an interpretation of the QAOA’s circuit as an interferometer operating globally in energy space. We will derive analytical results for the interference amplitude in later sections, focusing on the p = 1 single-layer ansatz (see Eq. 7): | Ψ QAOA - 1 γ , θ 〉 = R x θ ⊗ N e − i γ E ̂ H ⊗ N | 0 〉 ⊗ N . (10)

As sketched in Figure 1, the Hadamard gate splits the quantum state |0⟩ into a superposition of all states in the computational basis, acting like a multidimensional “mirror”: | s 1 〉 = | + 〉 ⊗ N = 1 2 N ∑ x | x 〉 . (11) Then, the evolution with the diagonal cost Hamiltonian E ̂ | x 〉 = E x | x 〉 , given by the operator U E ( E ̂ , γ ) , imparts phases on all states, acting like the branches of the interferometer: | s 2 〉 = e − i γ E ̂ | s 1 〉 = 1 2 N ∑ x e − i γ E x | x 〉 . (12) Finally, the mixing operator R _x (θ)^⊗N recombines the energy states so that the interference transforms the energy-dependent relative phases in measurable probability amplitudes: | Ψ QAOA - 1 γ , θ 〉 = R x θ ⊗ N | s 2 〉 = ∑ x F x E x , γ , θ | x 〉 . (13)

This intuition is much clearer when looking at the algorithm operating on a single qubit, with rotation R _x (θ): R x θ = e − i θ σ x = I ⁡ cos ⁡ θ − i σ x ⁡ sin ⁡ θ = cos ⁡ θ − i ⁡ sin ⁡ θ − i ⁡ sin ⁡ θ cos ⁡ θ . (14) Consider the following two-state subspace: | x 0 〉 = | 0 , x 1 , x 2 , … , x N 〉 , | x 1 〉 = | 1 , x 1 , x 2 , … , x N 〉 , where x _i = {0, 1} (so x ₀ and x ₁ differ by only 1 bit), with energies E x 0 ≡ 〈 x 0 | E ̂ | x 0 〉 and E x 1 ≡ 〈 x 1 | E ̂ | x 1 〉 . The probability amplitudes before the local interference are A x 0 2 = | ⟨ x 0 | s 2 ⟩ | 2 , A x 1 2 = | ⟨ x 1 | s 2 ⟩ | 2 , (15) respectively. After applying the local gate in Eq. 14 on the qubit that corresponds to the different bit, the interference shifts the amplitudes to N x 0 and N x 1 : N x 0 N x 1 = R x θ e − i γ E x 0 A x 0 e − i γ E x 1 A x 1 = e − i γ E x 0 A x 0 ⁡ cos ⁡ θ − i A x 1 e − i γ F x 1 , x 0 ⁡ sin ⁡ θ − i A x 0 ⁡ sin ⁡ θ + A x 1 e − i γ F x 1 , x 0 ⁡ cos ⁡ θ , (16) where F x 1 , x 0 ≡ E x 1 − E x 0 . Thus, the probability amplitudes become N 0 2 = A x 0 2 ⁡ cos 2 ⁡ θ + A x 1 2 ⁡ sin 2 ⁡ θ − A x 0 A x 1 ⁡ sin 2 θ sin γ F x 1 , x 0 N 1 2 = A x 0 2 ⁡ sin 2 ⁡ θ + A x 1 2 ⁡ cos 2 ⁡ θ + A x 0 A x 1 ⁡ sin 2 θ sin γ F x 1 , x 0 . (17) Note that the relative phase between the states γ F x 1 , x 0 = γ E x 1 − E x 0 controls the amplification. Therefore, assuming sin 2 θ < 0 and γ > 0 or vice versa and provided the angle γ is sufficiently small, the interference term always increases the population of the lowest energy state and reduces that of the other. Furthermore, Eq. 17 also reflects how the relative sign between sin 2 θ and γ controls the direction of optimization. The probability amplitude of the lowest/highest energy state is enhanced when the signs are opposite/equal.

This interferometric behavior expanded to the N-level scenario imposes non-trivial, structure-dependent upper bounds on the |θ| and |γ| angles to avoid symmetries and random scrambling of the energy states. In particular, Eq. 17 shows that for the single-qubit interference: n π ≤ θ < n + 1 π a n d | γ | < π F x max , x min ∼ O ‖ E ̂ ‖ − 1 , (18) with n ∈ Z . The N-qubit interference makes this picture richer and more complicated, as derived in the next sections. Analytical and numerical results for the single-layer QAOA on general Ising spin models show that the optimal |γ| for an N-qubit system actually scales as O ( N N | e | − 1 / 2 ) (Ozaeta et al., 2022; Díez-Valle et al., 2023), where |e| is the number of edges or non-null elements in the coupling matrix J (see Eq. 27). For instance, in the two-level system E ̂ = 1 2 Δ σ z , the probability of measuring the ground state |σ ^z = −1⟩ and the highest energy state |σ ^z = +1⟩ is | ⟨ σ z = ± 1 | Ψ γ , θ ⟩ | 2 = 1 2 1 + σ z sin 2 θ sin γ Δ , (19) which is maximal with optimal angles θ = ± π 4 and γ = π 2 Δ (see Figure 2).

Let us extend the two-level picture to a general framework dealing with the complete interference spectrum of a single-layer QAOA circuit. We focus on how the QAOA transforms the wavefunction amplitudes: F x ≡ ⟨ x | Ψ QAOA - 1 γ , θ ⟩ , (20) where x is any eigenstate in the energy spectrum that we will denote as the reference state. As previously shown, a local R _x rotation acting on the ith qubit mixes the amplitude between each state x with the corresponding state x′ that differs only on the value of the ith bit (see Eq. 17). Thus, the action of the whole mixing operator R x ⊗ N is that of a complete mixing between all states in the computational basis: F x = ∑ x ′ e − i γ E x ′ 2 N / 2 〈 x | ∏ j N cos ⁡ θ − i σ j x ⁡ sin ⁡ θ | x ′ 〉 . (21) The weight of each mixing process depends on the number of bits that must be flipped to go from one state x to the other x′, that is, the Hamming distance between both states, and also on the oscillating terms induced by the angles θ and γ. The amplitude of the state generated by the algorithm after the mixing operator U ( H ̂ x , θ ) is given by an interference formula, F x ≡ ⟨ x | Ψ ̃ ⟩ = 1 2 N / 2 ∑ x ′ cos θ N − H x , x ′ − i ⁡ sin θ H x , x ′ e − i γ E x ′ (22) where H x , x ′ is the Hamming distance between two bit configurations x and x′ that represent the eigenstates of a spin model with energy E _x and E x ′ . We can unify the weight sum terms into a single exponential by the following change of variables in which the rotation angle θ is reparameterized in terms of an exponent r and a normalization R, cos θ = R 1 2 ⁡ exp r 2 , sin θ = R 1 2 ⁡ exp − r 2 , (23) where, due to the symmetries of the interference, we consider θ ∈ ( 0 , π 2 ) without loss of generality. The sign of γ allows us to define if the interference increases the population of lower or higher energy states so that we cover all the possibilities. By this transformation, the interference amplitudes in Eq. 22 become a sum of exponentials over the entire configuration space, F x = R ⁡ exp r 2 N 2 ∑ x ′ exp − H x x ′ i π 2 + r − i γ E x ′ . (24)

Because Eq. 24 encompasses all eigenstates of the system encoded in its energy E x ′ and Hamming distance H x , x ′ with the reference state x, it is convenient to introduce a probability distribution relating distances in the computational base to the energy spectrum, p H , E ; x = 1 2 N ∑ x ′ δ H − H xx ′ δ E − E x ′ , (25)

This expression represents the relative number of eigenstates that, given a reference state x, have a Hamming distance H to that state and energy equal to E. In other words, this distribution captures the internal correlations inherent to the specific spin model. With this definition at hand, we can express the previous sum in Eq. 24 as the average over this probability distribution: F x ∝ exp r N 2 ∬ − ∞ ∞ ⁡ exp − H i π 2 + r − i γ E p H , E ; x d H d E . (26) Hence, the interference amplitude in Eq. 20 depends on the structure of the energy levels of the spin system manifested in the probability distribution p(H, E; x). In the following sections, we study such a structure to derive a unified expression of the QAOA probability amplitudes for universal Ising spin models.

Next, we focus on a universal set of non-deterministic polynomial time hard (NP-hard) Ising models (Barahona, 1982) and derive this scenario’s interference probability amplitude formula.

A broad spectrum of combinatorial optimization problems can be represented as a spin model described by an Ising Hamiltonian (Lucas, 2014), E I s = 1 2 ∑ i , j = 1 N s i J i j s j + ∑ i = 1 N h i s i , (27) where s = {−1,+1} ^N , N is the number of variables, J is an N-by-N square coupling matrix, and the magnetic field h is a vector of N coefficients. The quantum version of this Hamiltonian E ̂ I ( σ z ) is simply obtained by replacing the spin variables s with the corresponding Pauli-Z matrices σ ^z . The coupling matrix defines a structure in the problem that can be represented by a weighted graph with N vertices, connected by undirected edges i ↔ j that have associated weights J _ij , J _ji . In this work, we study families of models in which the J _ij coefficients are randomly drawn from a normal distribution N(μ = 0, σ ²).

This structure, together with the magnetic field values, defines families of optimization problems with different inner correlations and degeneracies. In the following subsections, we derive the interference amplitude for two scenarios that involve distinct energy level structures due to intrinsic global symmetries in the model. These differences necessitate a separate study of each case. Nevertheless, as explained below, the slightly different behaviors of these models converge to a common expression of the QAOA probability amplitude distribution |F(x)|² in Eq. 20. The two scenarios are represented by two well-known combinatorial optimization problems: quadratic unconstrained binary optimization (QUBO) (Kochenberger and Glover, 2004; Kochenberger and Hao, 2014), with a non-degenerate energy spectrum, and the maximum cut (MaxCut) (Nannicini, 2019; Sung et al., 2020; Harrigan et al., 2021), which exhibits a global Z 2 symmetry.