Band structure of two-dimensional graphene and of one-dimensional nanoribbons via variational quantum algorithms: from energies to wavefunction quality

Introduction

Strongly correlated electron materials represent a critical paradigm for future technologies due to the emergent physical phenomena driven by dominant electron-electron interactions[1]. Unlike conventional metals and semiconductors described by non-interacting band theory, strong Coulomb repulsion yields collective phase transitions, including metal-insulator transitions[2, 3, 4], unconventional superconductivity[5, 6], and spin-pairized ordering[1]. In carbon-based architectures, electron-electron interactions drastically renormalize electronic structures driving superconductivity in magic-angle and crystalline multilayer graphene systems[7, 8]. In reduced dimensions, such as graphene nanoribbons (GNRs), enhanced Coulomb repulsion and quantum confinement induce localized edge states along zigzag boundaries[9]. Furthermore, periodic structural engineering of GNRs gives rise to topological edge states and localized end states governed by symmetry-protected topological phases[10, 11]. Harnessing these interaction-driven superconducting phases alongside localized and topological boundary states require many-body accounting of electronic band structures to guide material engineering at the atomic scale.

Simulating strongly correlated many-body quantum physics on classical computers encounters severe computational bottlenecks. Exact diagonalization explicitly solves the $N$-body Schrödinger equation, but is strictly constrained to tiny clusters due to the exponential growth of Hilbert space dimension. Tensor network methods like Density Matrix Renormalization Group (DMRG)[12], deliver near-exact precision for quasi-one-dimensional systems by truncating low-entanglement states, however, extending DMRG to higher dimensional lattices requires an exponential expansion in bond dimension. Finally, Green’s function approximations and mean-field extensions such as $GW$[13] and Dynamical Mean-Field Theory[14] simplify electron correlations by mapping lattice interactions onto local impurity model, yet, they fail to capture non-local correlation fluctuations, long-range entanglement, and topological edge states without empirical parametrization.

Quantum computers circumvent classical simulation bottlenecks by using controllable quantum states to store and process physical wavefunctions[15]. By mapping fermionic operators onto qubit registers[16, 17], a system of $N$ spin-orbitals maps directly onto $N$ qubits, storing the full $2^N$-dimensional Hilbert space. Algorithms such as Quantum Phase Estimation (QPE)[18] leverage unitary time evolution to calculate ground-state energies and energy bands with polynomial circuit scaling. However, executing long-depth algorithms such as QPE remains infeasible on current devices. In the present NISQ[19] era, hardware is constrained by gate errors, crosstalk, and limited coherence times. To overcome these limitations hybrid quantum-classical approaches like the Variational Quantum Eigensolver (VQE)[20, 21] have emerged. By parametrizing low-depth quantum circuits and performing parameter optimization in classical hardware, variational algorithms minimize coherence requirements providing a practical route to calculate electron correlations and band structures[22, 23, 24, 25]. Furthermore, quantum strategies and scalable algorithms for correlated lattice systems such as the Fermi-Hubbard model[26, 27] and demonstrations of quantum utility[28] have expanded the reach of near-term processors in condensed matter physics.

Although the ultimate goal of quantum computing in materials science is to simulate strongly correlated many-body Hamiltonians beyond classical reach, implementing non-interacting tight-binding Hamiltonians on physical processors represents more than a preliminary step. Subspace algorithms, such as Quantum Krylov Subspace (QKS) methods, Krylov Quantum Diagonalization (KQD) [29], and Sample-Based Krylov Quantum Diagonalization (SKQD) [30], critically depend on high-quality initial reference states, as their convergence rates and subspace sampling efficiencies rely on the overlap with the true low-energy manifold. Likewise, warm-started Variational Quantum Eigensolver[31, 32] protocols leverage single-particle approximations to initialize variational parameters in favorable regions of the Hilbert space, mitigating barren plateaus and accelerating convergence toward correlated ground states. Single-particle tight-binding solutions offer a controlled starting point to build time-evolved Krylov bases or sample relevant electronic configurations or warm-start variational optimizers. In this context, Cerasoli et al. [33] pioneered mapping empirical tight-binding Hamiltonians of bulk silicon onto Pauli operators to solve for energy bands on IBM Quantum hardware, while Sherbert et al. and others [34, 35] established a generalized framework for calculating band structures from atomic orbitals. However, applying these algorithms to low-dimensional architectures such as one-dimensional graphene nanoribbons (GNRs) on physical quantum hardware remains largely unexplored.

In this work we address two carbon systems of different dimensionality within a single variational framework. On the one hand we take the two-dimensional monolayer graphene. On the other hand we take the quasi-one-dimensional GNRs. First, we perform the ideal simulation of the proposed method and we recover the exact diagonalization results. Second we perform simulations under quantum hardware noise models and we see that finite sampling and noisy optimization have comparable contributions. Third we perform executions in real IBM

Quantum processors and we demonstrate the feasibility of the computation on current NISQ devices. We also add a subspace containment analysis to address not only the energy, but also the wavefunction quality. The work is organized as follows: In section II we define the tight binding model of graphene and GNRs. In section III we define the quantum algorithm we will use to obtain the band structures of graphene and GNRs. In sections IV and V we show the band structure calculations. In section VI we show the results regarding subspace containment. Finally, we explain the conclusions in section VII.

Theoretical model

This section introduces the two systems studied in this work and the Hamiltonians from which their bands are obtained. We first describe two-dimensional monolayer graphene and then the quasi-one-dimensional nanoribbons, in both cases within the same tight-binding framework and in the Bloch representation, which reduces the problem, at each point of reciprocal space, to a finite-dimensional matrix $H(\mathbf{k})$. Throughout the paper, the exact reference against which every variational result is compared is the direct numerical diagonalization (exact diagonalization) of that Bloch matrix on a dense sampling of $\mathbf{k}$; in all band figures (Figure 3–Figure 9), the exact diagonalization result is drawn as continuous gray lines and the variational result as colored points.

Graphene is a hexagonal (honeycomb) lattice of carbon atoms with two triangular sublattices $A$ and $B$ (Figure 1a). In the tight-binding approximation, the $\pi$ electrons are described by

where $t$ is the first-nearest-neighbor hopping integral and $t'$ the second-nearest-neighbor one ($t' \approx t/28$). By periodicity, for each wave vector $\mathbf{k}$ the problem reduces to a $2 \times2$ Bloch Hamiltonian over the two sublattices,

with

where $\{\delta\}$ are the three first-nearest-neighbor vectors and $\{\delta'\}$ the six second-nearest-neighbor vectors. Their energies are $E_{\pm}(\mathbf{k})=c(\mathbf{k})\pm t\lvert f(\mathbf{k})\rvert$. At first nearest neighbors ($t'=0$) the diagonal term vanishes, $c=0$, and the bands are symmetric about the Fermi level, $E_{\pm}=\pm t\lvert f\rvert$, touching at the Dirac points $K$ where $f(\mathbf{K})=0$; the second-nearest-neighbor term $c(\mathbf{k})$ breaks that electron-hole symmetry and shifts the Dirac point in energy. The band structure is traced along the high-symmetry path $\Gamma\to K \to M \to\Gamma$ of the first Brillouin zone (Figure 1b).

We turn now to the one-dimensional systems. Confining the sheet yields a system that is periodic along a single direction (the ribbon axis), with a unit cell of $N$ atoms per sublattice. The Bloch Hamiltonian preserves the block structure

where now $f(k)$ is an $N \times N$ matrix and $k$ is the wave vector along the ribbon. The form of $f(k)$ depends on the type of edge [36] (Figure 2):

• Zigzag: main diagonal $-t(1+e^{ik})$; subdiagonal $-t$.

• Armchair: main diagonal $-t e^{ik/2}$; subdiagonal and superdiagonal $-t$.

The full Hamiltonian $H(k)$ has dimension $2N \times2N$; its $2N$ eigenvalues at each $k$ constitute the bands of the ribbon. In the implementation (Sec. III) we take $N = 2^{n_q-1}$, so that $H(k)$ has dimension $2^{n_q} \times2^{n_q}$ and is represented with exactly $n_q$ qubits.

Quantum methods

This section describes the computational workflow, common to both systems, that takes each Bloch Hamiltonian to a band structure obtained on a quantum processor. We first explain how the matrix $H(k)$ is rewritten as

an operator acting on qubits and which parametrized circuit prepares the trial state; we then describe the VQE and VQD algorithms used to obtain, respectively, the lowest band and the excited bands. We next define the three execution modes in which every calculation of this work can be run, discuss the computational cost of the noise simulation and the error channels it contains, and present the strategy of execution on real hardware. Finally, we introduce the subspace-containment diagnostic with which we validate the variational wavefunction.

For each $k$ the matrix $H(k)$ of dimension $2^{n_q} \times2^{n_q}$ is constructed and decomposed in the basis of Pauli-matrix products,

an operator acting on $n_q$ qubits. The trial state on which this operator is evaluated is prepared with a parametrized circuit: we use the hardware-efficient ansatz EfficientSU2 (with a number of repetitions that is stated for each calculation), on account of its low gate count and its suitability for NISQ devices.

The lowest-energy band at each $k$ is obtained with VQE, by minimizing the expectation value

by means of a classical optimizer. We have employed COBYLA, SPSA and SLSQP: SPSA is more robust against the stochastic noise of the measurement [37, 38] (preferable with finite shots and on real hardware), whereas SLSQP, being a gradient-based method, converges rapidly and accurately in ideal simulations.

The upper bands are obtained with VQD, by deflation: when computing the $j$-th state, a penalty term is added to the cost function that enforces orthogonality with respect to the states already found,

with a hyperparameter $\beta$, set to $\beta= 10$ throughout this work, large compared with the bandwidth $\sim6|t|$, so that the penalty dominates over any energetic gain from collapsing onto an already occupied level. To accelerate convergence along the path in $k$ we use the optimal parameters of one point $k$ as the seed for the next (warm-start). Two practical remarks about this deflation loop are worth recording. The first is that the sequential character of VQD makes the $j$-th state inherit the errors of the $j-1$ preceding ones and accumulate $j$ penalty terms, so that the quality of the reconstruction degrades systematically with the index of the state. The second is that, on a noisy cost function, a local optimizer started from a single seed can settle in a local minimum in which a state skips a level, it converges onto an eigenspace that is not its own, while the deflation penalty keeps it orthogonal to the previous states. We found that increasing the iteration budget does not correct this and can even aggravate it, since on a noisy cost function the optimizer keeps wandering around the wrong minimum, whereas restarting the optimization from a different seed does recover the correct eigenspace. We therefore optimize each state from $n_r$ seeds, the first being the warm-start one, the rest random, and retain the one of lowest cost, evaluated on the ideal simulator so that the choice does not depend on sampling noise. Unless stated otherwise, we use $n_r = 1$; $n_r = 3$ is used in the three-qubit containment calculations of § 1.5, where the fidelity of the excited states is itself the quantity under study.

Three execution modes are used throughout the paper:

• statevector: ideal simulation. The expectation value is evaluated exactly on the state vector, without sampling or noise.

• aer_noise: simulation with a noise model of a device (see below), with a finite number of shots.

• hardware: execution on a real IBM Quantum processor.

We do not add the ideal simulation with a finite number of shots. The reason is that the statistical error of a finite sampling is already contained in aer_noise, which adds the device errors on top of it, and, in the systems studied here, the results are indistinguishable from statevector within that statistical error. The decomposition of the error into its sampling and optimization contributions, reported below, is obtained instead from the two energies recorded in every mode: the ideal one (exact expectation value over the optimal parameters) and the measured one (with noise or hardware), which allows separating the robustness of the optimizer from the effect of noise on the energy.

For the aer_noise simulations we adopt the systematic rule of emulating a device with exactly the number of qubits required by the circuit with a fixed seed for reproducibility. This procedure can reduce the computational cost by a factor of 500 compared to the simulation using a simulated backend of 133 qubits and does not alter the noise model.

The emulated device has the native gate set $\{R_z, \sqrt{X}, X, CX\}$, and the noise model assigns to each of them a depolarizing channel with the corresponding gate error, composed with a thermal relaxation channel built from the gate duration and the coherence times $T_1$ and $T_2$ of the qubits involved; the measurement carries, in addition, a classical readout error (bit-flip on the recorded outcome). The virtual $R_z$ rotations are exact and carry no error, as on the real device. For the sizes used here the model yields one-qubit gate errors of $\sim10^{-4}$ (durations of 30–60 ns), two-qubit CX errors of $\sim1$–$4 \times10^{-3}$ (durations of 0.3–0.9 $\mu$s), readout errors of $\sim3$–$5 \times10^{-3}$, and coherence times $T_1,T_2 \approx120$–200 $\mu$s, values representative of the current generation of superconducting devices. Not included are correlated errors, crosstalk between neighboring qubits and the temporal drift of the calibration, all of which are present on the real machine. The model is therefore expected to underestimate the error with which a given circuit is measured on the device, mainly through the errors associated with two-qubit gates; § 1.4 quantifies this for the two-qubit circuits of the nanoribbons, while the single-qubit circuit of the monolayer, which contains no two-qubit gates, is essentially unaffected (Sec. IV).

Running the variational optimization directly on hardware consumes a lot of quantum computing time because the cost function must be evaluated at each step. We therefore adopt, in the hardware mode, a hybrid strategy: the ansatz parameters are optimized entirely in simulation and only the optimal circuit of each band and $k$ point is run on the real device to measure its energy. In this way the number of jobs is reduced to one per (band, $k$) pair: the hardware-measured energy of the optimal variational state. In all hardware calculations readout error mitigation is applied. This strategy, which we refer to below as the hybrid hardware strategy, is the one used in every real-device result reported in this work. The agreement in energy does not by itself guaranty that the state prepared by the circuit is the correct physical state. A variational state can approach the exact eigenvalue while carrying appreciable weight in neighboring eigenspaces. To validate the wavefunction we adopt the subspace containment diagnostic introduced in Ref. [39]. For a variational state $|\psi(\theta)\rangle$ and an exact-diagonalization eigenspace $\mathcal{H}_{\alpha}$, we define the subspace containment

where $\hat{P}_{\alpha}$ is the projector onto $\mathcal{H}_{\alpha}$. If $\{|\phi_{\alpha,i}^{\mathrm{ED}}\rangle\}$ is an orthonormal basis of that eigenspace, this becomes

Eq. (9) covers two cases. For a nondegenerate eigenstate the eigenspace is one-dimensional and the containment reduces to the usual squared overlap with the corresponding exact eigenvector: $C_{\alpha}=\left|\langle\phi_{\alpha}^{\mathrm{ED}}|\psi(\theta)\rangle\right|^2$, a number that says what fraction of the variational state is the physical state one intended to prepare, with $1-C_{\alpha}$ the weight leaked into the remaining eigenvectors. For a degenerate manifold, by contrast, the projection onto the full eigenspace is the appropriate quantity because it is basis-independent. In that case the $d$ variational states obtained sequentially by VQD within a $d$-fold degenerate manifold should be understood as $d$ states spanning the same eigenspace, not as energetically distinct levels, and what is meaningful is that each of them and, jointly, the subspace they span, resides within $\mathcal{H}_{\alpha}$.

In the systems studied here degeneracies are confined to isolated high-symmetry points: the Dirac point $K$ of the monolayer, where the two $\pi$ bands touch. The zone boundary $k=\pi$ of the zigzag ribbon, where the spectrum collapses onto three eigenspaces: $E=-|t|$ and $E=+|t|$, each $(N-1)$-fold degenerate, so that their degeneracy grows with the width, and the twofold eigenspace of the two edge states at $E=0$. Finally, the degeneracies of the armchair ribbon at $k=2\pi/3$ and $k=\pi$ (see Figure 7). At those points, the containment must be read over the manifold in the sense just described. Eigenvalues are grouped into eigenspaces numerically, with a tolerance $\mathrm{tol}_{\mathrm{deg}}=10^{-6}$ in units of $|t|$, so that the two regimes are separated automatically at each $k$.

The containment is evaluated on the noiseless statevector of the optimal parameters and on parameters obtained in the presence of noise.

Results: monolayer graphene (2D)

This section reports the band structure of monolayer graphene along the path $\Gamma\to K\to M\to\Gamma$. Being a single-qubit system, the monolayer is the simplest case and sets the reference against which the nanoribbons of Sec. V are read.

Fig. 3 compares the valence and conduction bands obtained with exact diagonalization and statevector mode, in the two approximations of the model. At first nearest neighbors (Figure 3a) the bands are symmetric about the Fermi level and touch at the Dirac point $K$; upon including the second nearest neighbors (Figure 3b) the electron-hole symmetry is broken and the Dirac point shifts in energy. In both cases the agreement with the exact reference is excellent along the entire path, which validates the VQE/VQD workflow for the two-dimensional system.

Fig. 4 compares the energy measured under noise with exact diagonalization. The optimizations are done with COBYLA and SPSA. In both cases the band structure is well reproduced, with mean deviations of $\bar{\epsilon}/t\approx5.8\times10^{-2}$ (COBYLA) and $\bar{\epsilon}/t\approx4.0\times10^{-2}$ (SPSA). The SPSA optimizer turns out to be more robust on average against the stochastic noise of the measurement, although both exhibit larger pointwise deviations in some regions of the path, in line with systematic comparisons showing that the relative performance of these optimizers is strongly problem dependent [40]. These deviations are consistent with the measurement noise level associated with a finite sampling of 1024 shots.

Fig. 5 compares the energy measured on the hardware with exact diagonalization along the same path. The two bands are reproduced with a mean deviation of $\bar{\epsilon}/t\approx2.1\times10^{-2}$ including the Dirac point at $K$, an excellent result that reflects that the circuit is shallow and the impact of device noise is small.

This error is smaller than the one obtained with the noise model in Fig. 4, which might seem to contradict the expectation that the model underestimates device noise (Sec. III). However, the two values measure different quantities. In a noise model, parameters are optimized under noise, meaning the total error $\bar{\epsilon}$ includes two distinct components: optimization error ($2.2\times10^{-2}$ for SPSA; $4.5\times10^{-2}$ for COBYLA) and final measurement error ($2.4\times10^{-2}$ for SPSA; $2.5\times10^{-2}$ for COBYLA). In contrast, the hybrid hardware strategy optimizes parameters in an ideal simulation, reaching exact energies within $10^{-7}t$, so its reported error consists solely of the final measurement component. When isolating pure measurement error, the physical device ($2.1\times10^{-2}$) and the noise model ($2.4$–$2.5\times10^{-2}$) agree closely. The device does not underperform the model here because this single-qubit circuit involves no two-qubit gates, eliminating unmodeled hardware noise like crosstalk and correlated errors, while single-qubit gate errors are negligible ($\sim10^{-4}$). Furthermore, readout error mitigation was applied on the real hardware but omitted in simulation. As a result, the remaining error in both cases is dominated simply by finite sampling noise.

Results: graphene nanoribbons (1D)

We turn now to the one-dimensional systems. This section reports, in order: (i) the validation of the method in statevector mode; (ii) the dependence of the spectrum on the ribbon width, obtained by exact diagonalization; (iii) the effect of the device noise model, including the decomposition of the error into its sampling and optimization contributions; (iv) the execution on real IBM Quantum processors; and (v) a scalability and error analysis as the number of qubits is increased, carried out in statevector and aer_noise modes.

Figure 6 compares the statevector results of GNBs with exact diagonalization. The band structure of both ribbons with $n_q = 3$ (width $N = 4$) over the entire interval $k \in[-\pi,\pi]$ is optimized with SLSQP using an EfficientSU2 ansatz with two repetitions. The agreement is excellent for the four bands in both ribbons. It is worth noting that reaching this precision in the excited bands required a sufficiently expressive ansatz (two repetitions) and a gradient-based optimizer (SLSQP); with a single-repetition ansatz the excited bands of the zigzag ribbon exhibit appreciable deviations.

We next examine how the spectrum depends on the width of the ribbon Figure 7 shows the exact band structure as the width is varied between $N = 4$ ($n_q = 3$) and $N = 64$ ($n_q = 7$). As $N$ increases, the number of subbands grows and the spectrum tends toward that of two-dimensional graphene. In the zigzag ribbons one can appreciate the nearly flat pair of edge-state bands at $E \approx0$ for $k \in[2\pi/3,\pi]$, which becomes flatter as the width increases. In the armchair ribbons a gap appears around $E = 0$ that varies with $N$. It should be noted that, by construction of the mapping ($N = 2^{n_q-1}$), the accessible widths are restricted to powers of two; a systematic study of the different width families of the armchair ribbons, with their well-known dependence of the gap on $N$ [41], would require a more flexible mapping.

To study the effect of noise we reproduce the band structure of both ribbons ($n_q = 2$, four bands) in the aer_noise mode with the SPSA optimizer. Figure 8 shows that the measured energy reproduces the exact band structure with good fidelity despite the noise. In units of $t$, the mean deviation with respect to the exact result is $\sim3.7 \times10^{-2}$ in the armchair ribbon and $\sim3.3 \times10^{-2}$ in the zigzag one (indicated in each panel of Figure 8). The decomposition of the error into its two contributions reveals that both are comparable: the measurement noise (measured energy versus the ideal energy of the optimal parameters) is $\sim2.5 \times10^{-2}$, consistent with the finite sampling of 1024 shots ($\sim1/\sqrt{1024}$ per term), while the effect of noise on the optimization (ideal energy of the optimal parameters obtained with noise versus the exact one) is $\sim1.5 \times10^{-2}$.

effect of noise on the optimization (ideal energy of the optimal parameters obtained with noise versus the exact one) is $\sim1.5 \times10^{-2}$. Overall, the SPSA optimizer maintains good precision despite the noise.

Next, we compute the complete band structure of the narrowest nanoribbons ($n_q = 2$, $N = 2$) on real IBM Quantum computers. We used the processors ibm_fez (armchair ribbon) and ibm_kingston (zigzag ribbon), both of the IBM Heron family. Figure 9 compares the energy measured on the hardware with exact diagonalization. In the armchair ribbon the four bands are reproduced with a mean deviation of $\bar{\epsilon}/t \approx4.4 \times10^{-2}$, distributed fairly evenly among the bands ($3.7$–$5.9 \times10^{-2}$). In the zigzag ribbon the mean deviation is $\bar{\epsilon}/t \approx6.5 \times10^{-2}$: the lowest and highest bands are reproduced with $2.1 \times10^{-2}$ and $2.3 \times10^{-2}$, while the two central bands concentrate most of the error ($1.4 \times10^{-1}$ and $7.2 \times10^{-2}$), with deviations of up to $0.29t$. These larger deviations are not related to the near degeneracy of the central bands at $E \approx0$: they are largest at small $k$, where the central bands are well separated from the rest of the spectrum, and smallest close to $k = \pi$, where they become degenerate.

Because the parameters are optimized in ideal simulation, where their energies agree with the exact ones to within $\sim10^{-4}t$, the whole hardware error is the error with which the device measures the optimal circuit. The same quantity in the aer_noise calculation of Figure 8, the measured energy compared with the ideal energy of the same parameters, is $2.6 \times10^{-2}$ (armchair) and $2.2 \times10^{-2}$ (zigzag), for circuits of the same size (two qubits, ansatz with one repetition, 1024 shots). The real devices therefore produce between 1.7 and 2.9 times the measurement error predicted by the noise model, unlike the single-qubit monolayer of Sec. IV. The hardware deviations are also systematic rather than statistical: in seven of the eight bands the measured energy is shifted on average toward the center of the spectrum ($E = 0$), the signature of incoherent errors that drive the prepared state toward the maximally mixed one, whereas in the noise simulation the mean shift of every band stays below $1.6 \times10^{-2}t$. This is the expected effect of the error sources that the model leaves out and that act mainly on two-qubit gates: crosstalk, correlated errors and the drift of the calibration (Sec. III).

Finally, we show a study of how the error and the resources scale as the ribbon width increases. We reconstruct by VQD the four lowest bands for $n_q = 2,\ldots,6$ (widths $N = 2,\ldots,32$), on a grid of seven $k$ points in $[0,\pi]$, with the SLSQP optimizer and an EfficientSU2 ansatz with two repetitions, held fixed across all sizes. The error $\bar{\epsilon}/t$ is the mean, over the four bands and the seven $k$ points, of the absolute deviation with respect to exact diagonalization.

The results show one favorable feature and one limiting feature. The favorable one is the scaling of the resources. The number of qubits grows logarithmically with the width, $n_q = \log_2(2N)$. The number of ansatz parameters is linear in $n_q$ ($6n_q$ for two repetitions); and thanks to the sparse character of the Bloch Hamiltonian, the number of Pauli terms grows approximately linearly with the width ($\sim N$ for armchair and $\sim2N$ for zigzag) rather than with the worst-case exponential $4^{n_q}$ of a dense operator (Figure 10b). Since the measurement cost of a variational algorithm is governed by the number of Pauli terms to be estimated, the quantum resources per $k$ point remain modest as the physical system grows.

$n_q$	$N$	Param.	Pauli $a$	Pauli $z$	Evals $a$	Evals $z$	$\bar{\epsilon}/t$ (statevector) $a$	$\bar{\epsilon}/t$ (statevector) $z$	$\bar{\epsilon}/t$ (aer_noise) $a$	$\bar{\epsilon}/t$ (aer_noise) $z$	$\bar{\epsilon}/t$ (hardware) $a$	$\bar{\epsilon}/t$ (hardware) $z$
2	2	12	3	4	2 566	3 617	$6.2 \times10^{-7}$	$1.6 \times10^{-7}$	$4.3 \times10^{-2}$	$3.1 \times10^{-2}$	$4.4 \times10^{-2}$	$6.5 \times10^{-2}$
3	4	18	5	8	7 315	11 427	$1.7 \times10^{-5}$	$4.1 \times10^{-4}$	$4.7 \times10^{-2}$	$8.3 \times10^{-2}$	–	–
4	8	24	9	16	12 207	30 150	$7.5 \times10^{-5}$	$1.8 \times10^{-2}$	$1.1 \times10^{-1}$	$1.0 \times10^{-1}$	–	–
5	16	30	17	32	15 823	60 768	$1.1 \times10^{-2}$	$1.7 \times10^{-2}$	$9.5 \times10^{-2}$	$5.7 \times10^{-2}$	–	–
6	32	36	33	64	19 170	121 074	$1.5 \times10^{-2}$	$7.0 \times10^{-3}$	$5.5 \times10^{-2}$	$6.1 \times10^{-2}$	–	–

The limiting feature is the precision at fixed ansatz depth. As Figure 10a shows on a logarithmic scale, the mean error grows by four orders of magnitude across the range studied, from $\sim10^{-7}$–$10^{-6}$ at $n_q = 2$ to $\sim10^{-2}$ for $n_q \gtrsim5$. The reason is that an ansatz of fixed expressivity has to cover a Hilbert space whose dimension doubles with every added qubit: with two repetitions the circuit has only $6n_q$ parameters, a number that grows linearly, and at some point the exact eigenvector ceases to be reachable.

Under noise the picture is different. The aer_noise error lies between $3.1 \times10^{-2}$ and $1.1 \times10^{-1}$ at every size and shows no systematic trend with $n_q$. Up to $n_q = 4$ it exceeds the statevector error by two to five orders of magnitude, so the intrinsic limitation of the ansatz is completely masked by the noise; only for $n_q \geq5$, where the ideal error has itself reached $\sim10^{-2}$, do the two come within an order of magnitude of each other. Separating the two contributions clarifies what happens. The optimization contribution (ideal energy of the parameters found under noise versus the exact one) varies between $2 \times10^{-3}$ and $8 \times10^{-2}$ without a clear trend. The measurement contribution (measured versus ideal energy of the same parameters) grows from $2.4$–$3.0 \times10^{-2}$ at $n_q = 2$–$3$ to $4.4$–$7.5 \times10^{-2}$ for $n_q \ge4$, and it also changes in nature. For $n_q \le3$ it is essentially statistical, with a mean shift of at most $1.2 \times10^{-2} t$, whereas for $n_q \ge4$ the measured energies lie systematically above the ideal ones, by $4$–$7 \times10^{-2} t$ on average, that is, shifted toward the center of the spectrum. This is the accumulation of gate errors in a circuit whose number of CX gates grows as $2(n_q - 1)$: the larger the circuit, the closer the prepared state gets to the maximally mixed one. The hardware columns, available for $n_q = 2$, must be read with the caveats stated in the caption. The hardware error ($4.4 \times10^{-2}$ armchair, $6.5 \times10^{-2}$ zigzag) contains only the measurement contribution, and yet it is of the same order as the total aer_noise error ($4.3 \times10^{-2}$ and $3.1 \times10^{-2}$), which contains both. Compared with the measurement contribution alone of the aer_noise calculation ($3.0 \times10^{-2}$ for both ribbons), the device is worse by a factor of 1.5 (armchair) and 2.2 (zigzag), in line with the analysis of Figure 9.

Recovering the precision of the ideal simulation requires increasing the number of repetitions, at the cost of a deeper circuit.

We repeated the reconstruction of the four lowest bands of the zigzag ribbon at $n_q = 4$ for one to four repetitions of the ansatz, in statevector and aer_noise modes, on the same grid of seven $k$ points and with the optimizer used throughout this work in each regime (SLSQP and SPSA, respectively); Figure 11 collects the result. In statevector the gain is dramatic: the mean error falls monotonically from $1.1 \times10^{-1}$ at one repetition to $2.6 \times10^{-5}$ at four, three and a half orders of magnitude, as the circuit becomes expressive enough to reach the exact eigenvectors. Under noise the trend is the opposite. The measured error is essentially unchanged between one and two repetitions ($1.00 \times10^{-1}$ and $1.04 \times10^{-1}$) and then grows, to $1.35 \times10^{-1}$ and $1.62 \times10^{-1}$ at three and four, while the two-qubit gate count of the ansatz goes from three to twelve. Separating the two contributions shows that the ideal energy of the parameters found under noise does improve with depth ($9.6 \times10^{-2}$ to $6.3 \times10^{-2}$), so the optimization is not what fails; what degrades is the energy actually measured, because a deeper circuit accumulates more gate error at every evaluation. At four repetitions the two regimes differ by nearly four orders of magnitude, and the expressivity that ideal simulation rewards is precisely what the device punishes. The optimal depth on a noisy processor is therefore the shallowest one the problem tolerates, one or two repetitions, which is what we use throughout, whereas in ideal simulation deeper is always better. Since our noise model underestimates the noise of the real device (§ 1.4), the turnover on hardware can only come earlier. The two error mechanisms in the NISQ regime, the limited expressivity of a shallow ansatz and the flattening of the optimization landscape by noise, are well documented [42, 43].

The two ribbons do not degrade alike. The zigzag ribbon loses precision earlier than the armchair one (it is already at $\bar{\epsilon}/t \approx1.8 \times10^{-2}$ at $n_q = 4$, where the armchair one is still at $7.5 \times10^{-5}$) and it demands more work from the optimizer: at $n_q = 6$ it requires 121 074 cost-function evaluations against 19 170 for the armchair ribbon, a factor of more than six. Both facts have the same origin: the lowest part of the zigzag spectrum is more crowded.

This crowding is not due to the edge states, which lie at $E \approx0$, in the middle of the spectrum, and for $n_q \ge4$ are not among the four lowest bands. It comes from the bulk bands: at the zone boundary $k = \pi$ the zigzag Bloch Hamiltonian has the eigenvalue $-|t|$ with degeneracy $N - 1$ (together with $+|t|$, also $(N - 1)$-fold, and the two edge states at $E = 0$), so the four lowest bands converge to the same energy as $k \to\pi$. On our grid, from $n_q = 4$ on, their mean spacing is less than half that of the armchair ribbon ($0.15 |t|$ against $0.37 |t|$ at $n_q = 4$, $0.011 |t|$ against $0.038 |t|$ at $n_q = 6$). In such a spectrum the deflation penalty of Eq. (7) has to separate states whose energies differ by very little, the cost landscape becomes correspondingly flat, and the optimizer needs more evaluations to resolve one level from the next.

The same crowding explains why the zigzag error decreases between $n_q = 4$ and $n_q = 6$ (from $1.8 \times10^{-2}$ to $7.0 \times10^{-3}$), a decrease that does not mean that the reconstruction improves. When a variational state is not the exact eigenvector, it mixes the target level with its neighbors or settles on one of them, and the resulting energy error is bounded by the spacing between those levels. As the ribbon widens, the four lowest levels get closer together: in the zigzag ribbon their mean spacing falls from $0.15 |t|$ at $n_q = 4$ to $0.043 |t|$ at $n_q = 5$ and $0.011 |t|$ at $n_q = 6$. A worse wavefunction can therefore produce a smaller energy error. Measured in units of the level spacing, the error grows monotonically with $n_q$ in both ribbons (zigzag: 0.12, 0.39 and 0.62 at $n_q = 4$, 5 and 6; armchair: 0.08 and 0.39 at $n_q = 5$ and 6). The energy alone thus underestimates how much the reconstruction degrades as the spectrum becomes denser, which is precisely the limitation that motivates the wavefunction-level diagnostic of the next section.

Wavefunction quality

The preceding sections show that the VQE/VQD workflow reproduces the energies of both systems with good precision. However, agreement in energy does not necessarily imply that the state prepared is the correct physical state. In this section we apply the subspace-containment diagnostic of Eqs. (8)–(9) to validate the reconstructed wavefunction. We proceed in the same order as before: first the monolayer, in two regimes statevector and aer_noise; then the nanoribbons, in a three-qubit calculation whose parameters are optimized under aer_noise; and we close with a joint discussion of the two systems.

For monolayer graphene, evaluating the subspace containment of the valence and conduction bands along the path $\Gamma\to K \to M \to\Gamma$ demonstrates high wavefunction quality in both noiseless and simulated noise regimes. As a single-qubit system, its reconstruction is exceptionally robust, maintaining $C_{\alpha} \gtrsim0.99$ for both bands across the Brillouin zone. At the degenerate Dirac point $K$, the manifold-projected containment approaches unity, confirming that the physical wavefunction is preserved in the two-dimensional system.

We turn now to the GNBs. In this case we have eight states instead of two, and three states to be reconstructed sequentially by VQD instead of one. Figure 12 presents the containment of the ground state and the first two excited states of the three-qubit ribbons, with the parameters optimized in the aer_noise mode. The reconstruction is good throughout, but not uniformly. The ground state is essentially exact in both ribbons ($C_{\alpha} \geq0.988$ over the entire path). The excited states remain high, but fall off with the index of the state and with the position along the path: in the zigzag ribbon they stay at $C_{\alpha} \geq0.91$, with the minimum (0.927 and 0.912) at $k \approx5\pi/6$, where the three lowest levels are already converging toward the threefold degenerate manifold at $E=-|t|$ that they form at $k=\pi$ and the spectrum becomes crowded; in the armchair ribbon they stay at $C_{\alpha} \geq0.97$ except for the second excited state at $k \approx5\pi/6$, which drops to $C_{\alpha}=0.809$, leaking 17% of its weight onto the eigenspace immediately above ($E=-1.031 |t|$, a mere $0.27 |t|$ away). This is the expected behavior of a deflation scheme: the last state of the sequence is the one that has accumulated the most penalty terms and inherits the imperfections of the two preceding ones, and it is therefore the first to leak toward its neighbor when the levels are close.

It should be stressed that these containment values are those of a converged optimization, in the sense of Sec. III: with a single seed per state, the same calculation produced spurious collapses ($C_{\alpha} \approx0.04$ at an isolated $k$ point of the armchair ribbon, $C_{\alpha} \approx0.55$ at one point of the zigzag one) which did not reflect any physics but a state that had skipped a level, and which disappear on optimizing from three seeds. We record this because it is a practical hazard of VQD on noisy cost functions and, in the absence of a wavefunction-level diagnostic, it would be easy to mistake for a genuine effect.

The degenerate points deserve separate reading. At the zone boundary ($k=\pi$) of the zigzag ribbon, the three lowest states form a threefold degenerate manifold at $E=-|t|$. This is a bulk manifold, not the edge states: for $N=4$ the spectrum at $k=\pi$ is $\{-|t|^{(3)},0^{(2)},+|t|^{(3)}\}$, so the edge states at $E=0$ are the fourth and fifth levels, which are not reconstructed here. The three variational states jointly reside in that eigenspace, with individual containments of 1.000, 0.998 and 0.959 over the manifold and a total close to unity. Almost all of the missing weight of the third state (3.8%) goes to the next eigenspace, the edge-state manifold at $E=0$. The armchair ribbon has its own degeneracies, a twofold eigenspace at $k=2\pi/3$, where the two excited states give 0.999 and 0.998 over the manifold, and a further one at $k=\pi$.

The contrast between the two systems is informative. In simulation, both are reconstructed correctly at the wavefunction level once the optimization has converged: the two-dimensional system stays at containment close to unity along the entire path, and the one-dimensional one keeps an essentially exact ground state and excited states above 0.9, with a residual leakage that grows with the index of the state and concentrates where consecutive levels approach one another. The diagnostic therefore does not, by itself, condemn the variational reconstruction under a device noise model; what it does is sort the difficulty, ground state before excited states, isolated levels before crowded ones, and expose optimization failures that the energy would have disguised.

Conclusions

We have computed, within a single VQE/VQD variational workflow and in the $\pi$-electron approximation, the band structure of two graphene systems of different dimensionality: the two $\pi$ bands of the two-dimensional monolayer and the complete set of subbands of the one-dimensional zigzag- and armchair-edge nanoribbons, from ideal simulation to real quantum hardware. In ideal simulation, the method reproduces the bands of both systems with excellent precision against exact diagonalization. Under noise, the band structure is recovered with deviations of a few hundredths of $t$, with comparable contributions from finite sampling and from the noisy optimization. Execution on real IBM Quantum processors demonstrates the feasibility of the computation on current NISQ devices. With the parameters optimized in simulation, the error with which the device measures the optimal circuit agrees with the noise model for the single-qubit monolayer, but exceeds it by a factor of 1.5–3 for the two-qubit circuits of the nanoribbons, where it also becomes systematic. The scalability study shows a favorable feature, the number of qubits grows logarithmically with the ribbon width, and a limiting one, the precision of a fixed-depth ansatz degrades as the dimension of the Hilbert space increases. Deepening the ansatz does not rescue that precision on a noisy device: a sweep over the number of repetitions improves the error by three and a half orders of magnitude in ideal simulation while making the energy measured under noise steadily worse, so that the useful depth is bounded from above by the hardware and not by the physics.

This work adds a validation of the wavefunction quality via subspace containment. Under a device noise model, the diagnostic shows that both systems are reconstructed correctly once the optimization has converged, containment close to unity for the monolayer, and above 0.9 for the three states of the nanoribbons, and that the residual degradation is ordered: it grows with the index of the state in the deflation sequence and concentrates where consecutive levels approach one another. The diagnostic also proved able to expose a failure mode of VQD on noisy cost functions, a state settling on an eigenspace that is not its own, which the energy alone disguises and which is corrected by optimizing from several seeds. Furthermore, it naturally captures the degenerate manifolds of both systems.

The workflow validated here transfers to a fully on-device optimization without modification. In a broader perspective, the methodology validated here on non-interacting electron systems constitutes a controlled starting point toward the study of topological and strongly correlated phases in realistic two and one-dimensional materials.

Acknowledgments

J.B. acknowledges financial support from the Department of Education, Universities and Research of the Basque Government (Eusko Jaurlaritza) and the University of the Basque Country (UPV/EHU) under Grant No. 2122-26. U. A. acknowledges support from the Gipuzkoa Quantum Programme financed by the Provincial Government of Gipuzkoa with project number 2025-QUAN-000021-01 and 2025-QUAN-000037-01. U. A. acknowledges support from the ELKARTEK programme financed by the Basque Goverment with project number KK-2025/00079. We acknowledge the use of IBM Quantum services for this work; the views expressed are those of the authors and do not reflect the official policy or position of IBM or the IBM Quantum team.

References