Engineering Reflective Metasurfaces with Ising Hamiltonian and Quantum Annealing

We present a novel and flexible method to optimize the phase response of reflective metasurfaces towards a desired scattering profile. The scattering power is expressed as a spin-chain Hamiltonian using the radar cross section formalism. For metasurfaces reflecting an oblique plane wave, an Ising Hamiltonian is obtained. Thereby, the problem of achieving the scattering profile is recast into finding the ground-state solution of the associated Ising Hamiltonian. To rapidly explore the configuration states, we encode the Ising coefficients with quantum annealing algorithms, taking advantage of the fact that the adiabatic evolution efficiently performs energy minimization in the Ising model. Finally, the optimization problem is solved on the D-Wave 2048-qubit quantum adiabatic optimizer machine for binary phase as well as quadriphase reflective metasurfaces. Even though the work is focused on the phase modulation of metasurfaces, we believe this approach paves the way to fast optimization of reconfigurable intelligent surfaces that are modulated in both amplitude and phase for multi-beam generation in and beyond 5G/6G mobile networks.


I. INTRODUCTION
T HE study of wave propagation mediated by metamaterials and metasurfaces (MSs) has been a longstanding topic in applied physics and modern engineering [1]- [8]. They offer a range of advantageous responses in beam steering [9], wavefront shaping [10], anomalous reflection/refraction [11], spatial processing [12] and computation [13], with practical applications in sensing, imaging, and wireless communication [14]. It has been shown that it is possible to control the reflection phase of metasurface groups of subwavelength elements accurately and flexibly to realize digital coding MSs. More recently, reconfigurable intelligent surfaces (RISs) have been widely studied in the wireless communication community [14]- [17]. While MSs are typically passive surfaces whose patterns are engineered to achieve a specific reflective/transmissive behavior, RISs are reconfigured via software to control the electromagnetic (EM) propagation dynamically, where the desired far-field wavefront can be achieved by tuning the local reflection phase of the elements.
Researchers have devoted substantial efforts in understanding how to devise pattern optimization methods to achieve a specific scattering profile, including genetic algorithms [18]- [20], impedance-based synthesis [9], electromagnetic inversion [21], statistical learning [22], as well as dynamical optimization via switching across profile states [23], [24]. Wavefront shaping has been demonstrated in the minimal case of binary reflection phase control, where a phase delay of either 0 or π can be impressed in the locally reflected field re-radiated by the metasurface element. For a square metasurface with N elements per dimension, a total number of N 2 degrees of freedom (DOFs) produce 2 N 2 possible phase configurations. A key question arises on how to efficiently select the phase configuration that matches the desired scattering profile. In other words, the enormous parameter space in metamaterial optimization needs to be explored quickly. Combinatorial exploration has been proposed in [25], [26], which flips a few elements at a time and test the improvement in the scattering profile to decide whether to keep or bring the phase back to the original state.
In this work, we propose to find the optimal phase configuration by a physics-based approach inspired by the quantum mechanical physics of spins. The Ising model is widely used in statistical mechanics to describe the spin state of arrays of quantum particles [27]. The formal analogy with the Ising Hamiltonian has been exploited fruitfully in protein folding [28], electromechanics [29] and photonics [30]. The latter case in particular associates the particle spin to the reflection phase of pixels in spatial light modulators (SLMs). Here, we dwell on this analogy and develop an Ising model for the metasurface array with prescribed scattering profile. The scattered wave energy is used to calculate the EM free-field Hamiltonian, which is found to configure as an Ising Hamiltonian. Thereby, the global solution of the problem, i.e., the values of the local reflection phases across the MSs/RISs, is obtained by computing the ground state of the equivalent Hamiltonian.
Another interesting aspect of this work is to leverage quantum computing (QC) to speed up the optimization of derived Ising Hamiltonian. In recent years, the remarkable progress made in QC hardware has defined a new, Noisy Intermediate-Scale Quantum (NISQ), QC era [31]- [33]. Here, we focus on quantum combinatorial optimization algorithms, which run on NISQ devices to search for an optimal solution over all the combinatorial states of MS/RIS elements. One well-known example is the Grover Adaptive Search algorithm [34] with provable quantum speedup, which has been recently applied to Constrained Polynomial Binary Optimization (CPBO) problems [35]. Another representative example is the quantum approximate optimization algorithm (QAOA) [36] using unitary operators and quantum superposition. It is noted that both above-mentioned algorithms run on universal gate-model quantum computers. There is another specialized analog computer, so-called quantum annealer (QA) [37], [38], which belongs to the adiabatic quantum computing (AQC) regime [39]. These quantum hardware received considerably interests lately due to the number of available qubits and programmability [40]- [44]. In this study, we explore the QA to efficiently search the ground state of Ising Hamiltonian model. The Ising coefficients are embedded into the physical grid of qubits on an advanced QA hardware, the D-Wave 2000-qubit (DW2Q) quantum adiabatic optimizer machine [45].
It is worth mentioning that there are classical algorithms that have been developed for Ising problems. In the literature, there are two types of classical algorithms used for the energy minimization of Ising models. One may recast the Ising model into mixed integer programming (MIP) model that can be solved by relaxation-based classical exact solvers. These include the branch-and-bound search [46] via semidefinite programming (SDP) or branch-and-cut search by incorporating cutting plane algorithm in the branch-and-bound scheme [47]. However, it is shown in [48] that such methods is less efficient than or at best comparable to quantum annealing. This is likely because of the need to solve the relaxation of the MIP model [49]. In recent years, many have shifted their interests to classical heuristic methods such as simulated annealing [50], [51]. Simulated annealing performs the energy minimization of Ising Hamiltonian via the implementation of Metropolis-Hasting algorithm on classical machines. There are several attempts to compare performance of these classical heuristic solvers to quantum annealing algorithms [52]- [54]. Yet, there is still no definitive advantage of choosing either categories as their performance depends on the problem size and complexity. In this paper, we shall limit our analysis to the QA. Rigorous numerical analysis and performance comparison between QA and classical heuristic algorithms will be considered for future work.
The rest of paper is organized as follows: In section II, we present relevant notations for the analysis of plane-wave scattering from a rectangular surface array. Secondly, we propose suitable Ising Hamiltonians that are used to represent the problem of achieving the desired scattering profile from the surface array. Based on the mathematical representation of Ising models, we introduce the quantum annealing algorithm to find the ground-state solution of Ising Hamiltonian, and discuss a generalization to larger-scale problem and higherorder modulated reflecting surfaces. In section III, we apply the proposed work on various problem sizes and prescribed scattering angles to evaluate the optimised solution and computational performance. Numerical results verify that the proposed framework promises a computationally efficient, general tool for the design and configuration of EM reflective surfaces. Finally, we draw some conclusions and future perspectives.

A. Problem Statement
Metasurfaces are artificial structures engineered with desired properties to control and manipulate EM waves. The specific structure considered in this work is a large planar metasurface array with configurable phase response in the element reflection coefficient. For the purpose of illustration, we start with the binary phase shift metasurface, in which the array elements exhibit 0 or π phase responses [3], [4], [19], [55], [56]. The extension to optimize higher order phase quantization will be discussed in Section II.E .
A generic problem statement of reflective metasurface in wireless application is illustrated in Fig. 1, including a transmitter (Tx), a receiver (Rx) and a passive reflecting metasurface array. Consider the case that the direct-link between Tx and Rx is blocked by an obstacle, the metasurface can be used to create a virtual line-of-sight link with enhanced signal reflection [57]. Assume a far-away Tx radiates a TE polarized uniform plane wave on the metasurface array with incident angle θ i . The size of the surface element is d × d and the total number of elements is M × N . Assuming a unit magnitude for brevity, the incident electric and magnetic fields can be expressed as: By using the physical optics (PO) approximation [58], we obtain the induced electric current on the element (m, n): where ψ (m, n) ∈ [0, π] is the phase of the reflection coefficient, and the reflection amplitude is assumed to be 1. There is a considerable literature on the design of such 1-bit binary elements. Interested readers are referred to [19], [59]- [61].
where k x =k sin θ s cos φ s It is thereby clear that the C (r) depends on the radial distance to the observation point, E s θ /E s φ and A s refer to the attributes of metasurface element factor and array factor, respectively.
Based on the scattered fields of Eqs. (4) and (5), the radar cross section (RCS) can be written as: 4πd 4 λ 2 |E s | 2 |A s | 2 (10) From Equation (10), we observe that the RCS is a product of single element scattering pattern |E s | 2 = |E s θ | 2 + |E s φ | 2 , and the power pattern due to the array factor |A s | 2 . Note that a rudimentary PO approximation of patch scattering is used for the sake of completeness in the RCS derivation. The goal of the study is to optimize the element-wise phase profile distribution, ψ (m, n), across the metasurface, such that the maximum scattered wave is directed towards the receiver.

B. Ising Hamiltonian Model
At this stage, we focus on the optimization of binary metasurfaces with 0 or π phase response. The array factor in Eq. (9) can be rewritten as: s mn e jkxmd e jkynd (11) where the discrete variable s mn = ±1 corresponding to the 0/π element phase modulation, and the subscript indicates the element index (m, n). The scattering intensity of the array factor is obtained as: The expression in Eq. (12) naturally reminds us of the classical Ising spin lattice model, which has been originally introduced as a mathematical model for understanding ferromagnetism in statistical mechanics. The Ising model consists of discrete variables that describe the magnetic moments of atomic spins. Each spin can take one of two possible states: +1 for spin up and -1 for spin down. In the Ising model, the energy of the system is conveniently expressed as an effective Hamiltonian function. In what follows, we will discuss how to construct the target EM energy Hamiltonian for achieving the desired scattering profile from the metasurface array.
1) Desired Signal Maximization: One typical application of reflecting metasurfaces in wireless communication is the directional signal enhancement at the Rx. Namely, a judicious selection of element coefficients s * * leads to a constructive beamforming in the desired direction towards the Rx, thus improving the data transmission performance. The corresponding mathematical representation is: where the power patterns of element factor |E s | 2 and array factor |A s | 2 are functions of (θ s , φ s , θ i ). The solid angle integration in Eq. (13) is defined as: in which the desired direction is denoted by φ r and θ r , and ∆ φ and ∆ θ are introduced to account for the finite beamwidth. We will utilize the 3-dB beamwidth derived in [62]. The maximum optimization in Eq. (13) is equivalent to finding the groundstate solution of the effective Hamiltonian: (15) In the analogy with Eq. (12), the variable s * * takes values of +1 or −1, which can be interpreted as the spin value in the Ising model. The spin-spin interaction strength, J r mniv , is defined as: with C > 0 a positive constant. Here we choose C = 1/S r to normalize the energy Hamiltonian. Note that we have included a minus sign in Eq. (16) such that the lowest energy state of Ising Hamiltonian corresponds to the maximum scattered power. The extension to the case of multiple Rx destinations is rather straightforward by incorporating multiple solid angle integrations in Eq. (14).
2) Interference Suppression: Often in practice, wireless communication scenarios require the scattered wave to avoid certain areas or angular domains, in order to suppress the cochannel interference from nearby base stations, or to protect safety critical electronic devices. These attributes can be recast into an argmin optimization problem: The solid angle integration in Eq. (17) is defined as: specifies the angular domain where the destructive beamforming is desired.
The optimization in Eq. (17) can be recast into finding the ground-state solution of the effective Hamiltonian: The spin-spin interaction strength, J c mniv , is defined as: In the most general case, we may combine both effects of desired signal maximization and interference suppression, thus leading to a preferred signal hotspot as well as interferencefree zone [57]. This can be solved with an energy functional consisting of two components: where H r and H c are defined in Eqs. (15) and (19), respectively. The ground-state solution (lowest energy) of the Hamiltonian (21) will satisfy both requirements. It is noted that the computational complexity for calculating Eqs. (16) and (20)  C. Quantum Annealing Optimization 1) Introduction to Quantum Annealing: So far, we have recast the problem of achieving the desired scattering profile into finding the ground-state solution of Ising Hamiltonian for a M × N spin system. The spin-like degrees of freedom are offered by the discrete phase values achievable by the tunable metasurface element. Clearly, due to the enormous design space, 2 M N for M × N Ising spins, finding the optimal solution with classical computational algorithms can be very challenging. This is where the quantum annealing becomes an appealing way forward.
The essential concept of quantum annealing (QA) is to harness the natural evolution of quantum states to solve the energy minimization problem represented by the Ising spin glass model [63], or equivalently, the quadratic unconstrained binary optimization (QUBO) problem [64]. The QA starts from the quantum-mechanical superposition of all possible states (candidate states). Then the system evolves according to the time-dependent Schrodinger equation. Physically, this is realized by introducing a controllable quantum transverse-field term into the Hamiltonian that decreases adiabatically. Thus, the range of sampled solutions shrinks as the field strength decreases over time until the low energy state Hamiltonian remains. Interestingly, the QA can escape local minima solution via cooperative quantum tunneling effect. This is comparable to similar features of classical simulated (thermal) annealing (SA) [65] to escape local minima by thermal activation. Hence, the quantum annealing is a global heuristic optimizer, like simulated annealing, to find the global minima solution.
2) D-Wave Quantum Annealing Process: To achieve the ground state efficiently, we have compiled the Ising models in Eqs. (15) and (19) into a physical QA hardware, the D-Wave 2000Q (DW2Q) quantum annealer device. The DW2Q has 2048 functional quantum bits (qubits) represented by circulating currents in superconducting loops.
To demonstrate the QA algorithm, we consider an Melement linear metasurface array with the application of desired signal maximization. The corresponding quantum Hamiltonian with Ising spins in a transverse field is given by [66]: (22) whereσ x,z m are the Pauli spin matrices. Theσ z m represent the spin projections along either the +z or −z direction (taking values +1 for spin up and -1 for spin down). The A(t) represents the transverse Hamiltonian due to an applied transverse field in the x-direction.
The quantum annealing process starts at time t = 0 with A(0) B(0). The ground state of the initial Hamiltonian is given by the product state of the spins in the x-direction, which is an easy state to set up and initialize. The system is then evolved by decreasing A(t) and increasing B(t) until the annealing time t f is reached. If the increasing in B is slowly enough, by the adiabatic theorem the final state of the system is the ground state of the target Hamiltonian. Namely, the qubits have dephased to classical systems, and theσ z m can be replaced by classical spin variables,ŝ 1 , · · · ,ŝ M , which represent the ground state (optimal) solution.
3) Embedding Ising into DW2Q Hardware: It is worth noting that the qubits of DW2Q quantum processing unit (QPU) are not fully connected. Instead, the interconnection of qubits is represented by a topology known as the Chimera graph. The Chimera graph comprises of a 16 × 16 two-dimensional lattice of unit cells. Each unit cell consists of four horizontal qubits and four vertical qubits. The couplings among qubits within a unit cell are denoted as internal couplers, while the couplings between different unit cells are referred to as external couplers. A cropped view of 2 × 2 grid of unit cells in Chimera graph is illustrated in Fig. 2, from which we can see that any given qubit is connected with at most six other qubits.
Recognizing that the Ising variables in Eqs. (15) and (19) are almost all-to-all connected, thus, a process of mapping the fully connected graph of our Ising problem into the Chimera graph is required. This is achieved by embedding each Ising variable (hereafter denoted as logic node) onto a connected chain of physical qubits presented in the Chimera graph [67]. In other words, the chain is a collection of physical qubits bound together to represent a single logical node. The physical qubits within a chain are constrained to have the same value, thus acting like a single logical node.
Consider a metasurface with 8 Ising variables as an example. It gives rise to a fully connected (complete) graph with 8 logic nodes, as shown in Fig. 3(a). A uniform triangle embedding [68] of this complete graph to the Chimera graph is depicted in Fig. 3(b). The color is introduced to reflect the mapping from logical node to physical qubits. Clearly, each logic node in the complete graph is represented by a chain of 3 physical qubits in the Chimera graph. The details of the triangle embedding are discussed in [68] for interested readers.  Note that due to the constrained connection in Chimera graph, the chain length will increase accordingly with the size of the fully connected graph. For the case of N Ising variables, the length of the chain needs to be (N/4 + 1), thereby a total number of N 2 /4 + N physical qubits are required in the embedding process. Since the maximum available physical qubits in DW2Q is 2048, a straightforward calculation shows that the largest size of fully connected graph can be solved by the Chimera architecture is 65 using the triangle embedding.
D. Further Technical Discussion 1) Converting into Non-zero Bias Ising Model: As we compare Eqs. (15) and (19) to typical Ising problems in statistical mechanics, they appear to have only quadratic coupling coefficients, denoted by the spin-spin interaction strength J . The linear coefficients (so-called biases) for Ising spins are all zero. Such zero-bias problems are known to be harder to solve since there are many degenerate states [69]. As an example, for any ground-state solution, we would obtain the same energy state by reversing the spin value of all Ising variables.
We present an easy way to convert Eqs. (15) and (19) into non-zero bias Ising models. For the purpose of elucidation, we use an M -element linear metasurface array along thex-axis as an example. The Hamiltonian in Eq. (15) will reduce to: The first step is to fix the spin value of a random Ising variable. Without loss of generality, we set the 1 st Ising variable to be 1. Then the Hamiltonian in Eq. (23) can be rearranged as: where J r 1(m+1) are the biases and J r (m+1)(i+1) are the couplings between spins. Equation (24) thereby becomes a standard non-zero-bias Ising problem, whose ground state solution is same as Eq. (23).
2) Two-Level Optimization for Large Array: Due to the mismatch of the fully connected graph generated from our Ising models and the very sparse Chimera graph in D-Wave QPU hardware, the size of Ising problems can be solved is relatively small (around 65 Ising variables). On the other hand, a realistic 2D metasurface may have hundreds or thousands of Ising variables. To overcome this QPU hardware limitation, we will discuss a two-level QA optimization for larger-scale metasurface problems.
The idea is inspired from the subarray concept in modern phased array systems, in which antenna elements are grouped into subarrays and then subarrays form the entire array. Consider the M by N metasurface shown in Fig. 1 We first write the Hamiltonian of a M (1) by N (1) subarray configuration: where p mn and p iv are element spin variables, and the spinspin interaction is a weighted integral by the angle-dependent, element scattering intensity, same as the one given in Eq. (16). By utilizing the QA algorithm described in the previous subsection, an optimized sequence of Ising spin variables, p 11 , · · · ,p M (1) N (1) , can be obtained. Subsequently, we can write the optimized subarray scattering intensity |Ẽ s | 2 as: The next task is to superimpose the M (2) by N (2) subarrays to achieve a desired signal maximization. Note that M (2) = M/M (1) and N (2) = N/N (1) . Each subarray can take one of two possible states: spin up (+1) or spin down (-1). Accordingly, we can write the Hamiltonian for those subarray contributions as: where the subarray spin-spin interaction is calculated as: (28) where the D x and D y are subarray center-to-center spacing along thex-axis andŷ-axis as shown in Fig. 4. The computational complexity for calculating the interaction in (28) scales with the number of subarrays.
After the optimized sequence of subarray spin variables, q 11 , · · · ,q M (2) N (2) , is obtained, the final Ising variable for individual array element can be determined byŝ mn =p mn ·q mn in a post-processing step. The computational complexity in this post-processing step scales with the number of array elements.
It is noted that the concept can be compared to the renormalisation group approach in computational quantum physics [70], [71]. The subarray can be viewed as block spin in the scale transformation.
3) Extension to Quadriphase Modulation: So far, the discussion of this paper is focused on the binary modulated metasurface array, where individual elements have 0 or π phase responses in the reflection. Such binary elements have been widely studied in the literature due to the low hardware complexity and cost. Nevertheless, the low resolution in the phase response also restricts the maximum achievable directivity of RCS pattern [72]. In this subsection, we will extend the proposed work to the quadriphase modulation.
Recall that the construction of Ising Hamiltonian starts by representing the reflection coefficients with Ising spin variables. Consider the metasurface element with four phase responses, ψ (m, n) ∈ {0, π/2, π, 3π/2}. The complex reflection coefficients Γ (m,n) = e jψ(m,n) are given in Table 1. ψ (m, n) 0 π/2 π 3π/2 e jψ(m,n) 1 + 0j 0 + j −1 + 0j 0 − j It is clear that both real and imaginary parts of Γ (m,n) have three states {−1, 0, 1}. The direct implementation requires a linear combination of four Ising spin variables. Nevertheless, we can rewrite the reflection coefficient as: Γ (m,n) = e jψ(m,n) e −j π 4 = s Re mn + js Im mn e −j π 4 √ 2 2 (29) where only two Ising spin variables s Re mn ∈ ±1 and s Im mn ± 1 are used to represent the four phase states in the reflection coefficient. The detailed mapping is given in Table II. ψ (m, n) π/4 3π/4 5π/4 7π/4 e jψ(m,n) By substituting Equation (29) into Eq. (9), we obtain the scattering intensity of the array factor as: As indicated by Eq. (30), the power pattern of array factor is now expressed as a quadratic polynomial of Ising spin variables. In the case of desired signal maximization in Eq. (13), we can transform the argmax into the equivalent Ising form with the Hamiltonian defined by: The spin-spin interaction strengths are determined by: Compared to the binary phase configuration, the total number of Ising variables is increased by a factor of 2, namely, 2M ·N for a M by N metasurface array. After the optimized value of Ising spin variablesŝ Re mn andŝ Im mn are obtained, the metasurface element phase response can be retrieved by: IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, 2021 7

III. NUMERICAL EXPERIMENTS
The main goal of this section is twofold: to validate the effectiveness of Ising Hamiltonian for achieving desired scattering profiles, and to study the performance of quantum annealing to find global optimization in a large search space. The section starts with the optimization of binary phase shift metasurfaces for linear and rectangular arrays. Subsequently, we present results for quadriphase metasurfaces, after which, the performance of quantum annealing is evaluated with regard to generic success metrics and problem scaling. The validation with the full-wave simulation of a square patch metasurface is presented to conclude the section.

A. Binary Linear Array Configuration
The first experiment is performed for a N -element linear metasurface array. The size of the array element is d = 1λ. The N elements are placed along theŷ-axis. The TE polarized plane wave is normally incident upon the array.
We would like to optimize the phase profile of the metasurface array, such that the scattering power towards a desired direction θ r is maximized. Concerning a linear array of size N = 8, the QA optimized Ising spin values for various preferred angles, θ r = 10 o , 20 o , 30 o are presented in Table III. The results using the brute force search are also supplied for comparison. As is evident, the QA optimization successfully found all optimal solutions corresponding to the ground-state solution of the Hamiltonian. Next, we plot the RCS pattern as a function of the scattered angle θ s in Fig. 5. We can see that, by using the optimized Ising spin values as the metasurface phase profile, the main beam of the RCS pattern indeed is directed towards the preferred direction. The results validate the effectiveness of the proposed Ising Hamiltonian for desired signal maximization.  Next, we keep the same metasurface element setting, and vary the array size by N = 6, 8, 10. The optimized Ising spin values in comparison with the brute-force search results are shown in Table IV. Again, the same optimal solutions are found. The resulting RCS patterns are depicted in Fig. 6. We observe that the increased array size leads to a narrower, more directional main beam towards the desired direction θ r = 20 o .

B. Binary Rectangular Array Configuration
The second example considered here is a 4 by 4 rectangular array with binary phase response in the element reflection coefficient. We choose a relatively small array in the study such that a bruce-force search approach is afforded for comparison. Same as the linear array case, we set the size of the array element as d = 1λ, and the TE polarized plane wave is normally incident upon the array. Figures 7 and 8 show the RCS pattern and optimized Ising spin value for the desired direction θ r = 20 o , φ r = 40 o , and θ r = 30 o , φ r = 150 o , respectively. In both cases, the scattered power is maximized at the desired direction. The optimal sequence of Ising spin values are provided in Table  V, in comparison to the bruce-force search results. The results verify the effectiveness of our Ising Hamiltonian model and QA optimization algorithm in 2D planar array configurations.   , we can include the argmin Hamiltonian using Eqs. (19) and (21). In Fig. 9(a), we plot the RCS heatmap using Ising Hamiltonian for desired signal maximization at θ r = 30 o , φ r = 150 o . The result after incorporating the interference suppression in the energy Hamiltonian is shown in Figure 9(b). We clearly observe that the scattered power is significantly smaller in the prescribed interference-free angular domain.

C. Quadriphase Array Configuration
We now consider reflective metasurface arrays with quadriphase element response. The element reflection coefficient has four quantized phase states with a π/2 phase increment, therefore it is often referred to as 2-bit reflecting element. The construction of the quadriphase Ising Hamiltonian has been discussed in details in Section II. D(3). Note that 2 Ising spin variables per element are required in this case, and the element phase response can be retrieved using Eq. (36). Figure 10 shows the RCS pattern of a 1D 2-bit metasurface array for desired signal maximization at θ r = 20 o , where the results of the binary (1 bit) metasurface array are also included for a quantitative comparison. We notice that the quadriphase configuration significantly increases the maximum achievable scattering power at the desired direction.
Next, we perform the experiment on a 4 by 4 quadriphase metasurface array for the desired direction θ r = 20 o , φ r = 40 o . The RCS pattern and optimized phase profile are shown in Fig. 11(a). Comparing to the result of binary array in Fig.  7, more than two times more scattered power is observed with the quadriphase array. Finally, we keep the same array setting, but change the incident angle from normal incidence θ i = 0 to oblique incidence θ i = 30 o . In a wireless communication setting, this corresponds to the case that the Tx has moved to a new location while the Rx is stationary. The results are shown in Fig. 11(b), as it is clear the main beam of the scattered power is directed towards the desired direction.

D. Performance Study of Quantum Annealing
As compared to universal gate-model quantum computers, the D-Wave quantum annealer hardware belongs to the adiabatic quantum computing (AQC) regime, which is specialized in solving NP-hard combinatorial optimization problems. Due to the undesired noise on the quantum processors, it is difficult to fulfill the precise adiabaticity condition in practice. Therefore, there is no theoretical guarantee to find the ground state (i.e. an optimal solution) for every single run. But the global optimality may be reached with high probability.
Given the probabilistic nature of quantum annealing, we will evaluate the performance with respect to two metric: 1) the probability of producing an optimal solution in a singleinstance run, and 2) the time-to-solution (TTS) required to find the ground state at least once with some desired probability.
1) Metric in terms of Success Probability: To evaluate the single-run success probability of the QA algorithm, we repeatedly run the annealing procedure many times, collect the return energy sample from each independent run, and then plot the discrete probability distribution of the energy samples. The default parameters in the D-Wave quantum annealer [73] (annealing time = 20µs, chain strength=1.0, heuristic minor-   Figure 12 shows the probability distribution for the energy samples collected from 400 runs of a N =8 linear metasurface array problem. The ground state energy is referred to as E g and the energy sample is denoted by E. The probability of reaching the ground state (global optimal solution) is around 95%, which is often referred to as the success probability of a single-instance run, P s [53], [54]. We then perform the experiment for the 4 by 4 binary metasurface array with energy Hamiltonian for desired signal maximization as in Fig. 7. The result of this problem instance, denoted as problem A, is presented in Fig. 13. We note that the probability of success of reaching the ground state is almost 1. Next, we run the experiment for the same binary metasurface but with energy Hamiltonian H r + H c as in Fig. 9(b). The probability distribution of this problem instance, denoted as problem B, is depicted in Fig. 14. The success probability of finding the global optimal solution is now around 75%.

Optimal solution
To investigate the problem-dependent success probability, we repeat the two problem instances with the brute-force search approach. After all samples are collected, we sort the solutions with respect to the normalized energy E/E g . The first 50 solutions for problems A and B are plotted in Fig. 15(a) and (b). Generally speaking, the QA success probability is related to the energy landscape of the Ising Hamiltonian. If the energy gap between near-ground states and ground state is small, e.g. Problem B in Fig. 15(b), there is a moderate possibility that the QA algorithm returns a nearoptimal solution, as shown in Fig. 14. Whereas for the case of distinguishable low energy states, the QA can find the ground state solution with very high probability in Fig. 13.
In addition, we notice that for RIS applications an approximate global optimum is sometimes desired instead of a precise global optimality. In this case the near-ground states from QA runs can also be of some value.
2) Metric in terms of Time-to-Solution: In this subsection, we are interested in studying the performance of QA algorithm over a range of problem sizes (# Ising variables). As is well known, the combinatorial optimization task becomes increasingly complex as the number of Ising variables grows, since each additional Ising variable doubles the number of states over which the energy landscape is defined.  Fig. 9(b)).
The linear binary metasurface array is used as the test problem. To account for the problem-dependent energy landscape, we create three problem instances corresponding to θ r =10 o , 20 o , 30 o , collect energy samples from 400 runs of each problem instance, and then calculate the average success probability. The results for varied array sizes from N = 10 to N = 60 are depicted in Fig. 16. The success probability of finding the ground-state solution in a single run, P s (N ), decreases sub-linearly as the problem size increase. Note that if we are interested in finding the ground state at least once with some desired probability P d , we can calculated the required number of runs by: Taking the problem N =60 as an example, with 17 independent runs we can achieve 99.99% probability of finding the ground state solution. The other results are shown in Fig. 16. We can then define the time-to-solution metric to achieve a given target probability P d , which requires R(N ) number of QA runs as introduced in Eq. (37). During our access to the D-Wave quantum processor, we record the QPU access time including QPU initialization, programming, and QA sampling time (anneal/readout/delay time). As is seen in Fig. 17, the QPU access time grows linearly with increasing problem size, which is an encouraging result using the QA optimization for metasurfaces. In addition, we also plot the CPU execution time that is needed to compute the spin-spin interaction terms in Eq. (16), which also scales linearly with respect to the array size, i.e. O(N ).

E. Full-wave Simulation as Validation
We conclude the section of numerical experiments with the full wave simulation of a square patch metasurface. As shown in Fig. 18, each metasurface element consists of 7 by 7 metallic patches printed on a dielectric substrate above a metallic ground plane. The thickness of the substrate is h = 2mm with relativity permittivity of 4.0. The size of the sub-wavelength unit cell a = 5mm. Thereby the metasurface element has the size of d = 35mm, which is equal to one wavelength at the operating frequency 8.57GHz. By controlling the size of the metallic patch, p, we can achieve the desired quadriphase shift in the element reflection coefficient, as illustrated in Table VI.   The experiment is performed on a 16 by 16 quadriphase metasurface array with the normal incident plane wave and the desired direction θ r = 20 o , φ r = 40 o . Since the number of Ising variables, 512, exceeds largest fully connected graph solvable in the Chimera hardware, we have used the two-level optimization approach detailed in Sec. II.D. The problem is decomposed into 4 by 4 subarrays and each subarray includes 4 by 4 array elements. The QA optimized phase profile and RCS pattern are shown in Fig. 19(a) and 19(b). It is clearly seen that the scattered power is maximized at the desired direction. The power values of the main lobe peak and highest side lobe peak are 56.19 dBsw (at θ s = 19.9 o φ s = 40 o ) and 47.47 dBsw (at θ s = 16.1 o φ s = 52.3 o ), respectively.
We proceed with the full-wave simulation of the same metasurface array using the geometry-aware domain decomposition method [74], [75]. The computed RCS is depicted in Fig.  19(c), where a similar pattern is observed comparing to Fig.  19(b). The power values of the main lobe peak and highest side lobe peak are 55.42 dBsw (at θ s = 20 o φ s = 40 o ) and 46.85 dBsw (at θ s = 16 o φ s = 52 o ), respectively. Note that the absolute scattering power of the main beam is a bit smaller due to the non-ideal element-scattering amplitude. Still, the results validate the applicability and potential of proposed work for practical reflective metasurfaces.

IV. CONCLUSION
Over the last few years, we have witnessed an extensive and growing interest in leveraging reconfigurable intelligent surfaces towards smart wireless environments. One key question arises on how efficiently to select the phase configuration that produces a scattered field matching the desired scattering profile. This is of paramount importance when a solution to the optimization problem is not available in closed-form, and thus constitutes a substantial computational task.
Whereas recent researches focus on artificial intelligence and deep learning, this paper takes another direction aiming to leverage the power of quantum computers to overcome the computational optimization complexity. The scattered wave power is expressed as an Ising Hamiltonian: a common mathematical abstraction employed in statistical mechanics to describe the spin state of arrays of quantum particles. An analogy can be made between the discrete meta-atom state and the spin degree of freedom in order to design the reflection phase mask of metasurfaces and RISs.
The results show a viable way forward for analyzing and controlling the interaction of large reconfigurable surfaces and complex radio environments. This work constitutes the first step towards including quantum computation within EM optimization tools for the next generation of wireless mobile networks. Future work includes multi-level spin optimization for multi-phase RISs, joint optimization of reflection amplitude and phase modulation as well as the incorporation of mutual coupling effects between array elements. The colormap from black, dark gray, light gray to white corresponds to the element phase state ψ (m, n) = 0, π/2, π, and 3π/2.