The model on this page is the full, spinful two-dimensional Fermi–Hubbard model with hopping t, diagonal hopping t′ and local repulsion U. Its parameters are U/t = 8 and t′/t = −0.25. Here “full” means retaining every term of this specified model and the fermionic signs. It does not include every possible interaction in a real material. The historical bare-XY pilots and the later full fermionic route are different calculations; results from the former do not validate the latter.
The explanation below starts with lattice sites and occupation. It then develops the energy and time units, the qubit representation and an exactly solvable example. Throughout, we distinguish the physical model, its representation on qubits and the accuracy of an execution.
1. What do 2D, XY and the dimension of the state space mean?
2D describes the geometry: each lattice site has two coordinates, x and y. Our square 6 × 6 lattice has 36 sites. Electrons can move horizontally and vertically; t′ adds diagonal connections. Open boundaries mean that opposite edges are not connected.
i=x+L_x y,\qquad 0\leq x\lt L_x,\quad 0\leq y\lt L_y,\qquad L=L_xL_y=36.
\]
XY describes operators: an XY coupling contains Pauli X and Y operators, for example XX + YY. These letters are not lattice coordinates. An XY spin model can itself live on a two-dimensional lattice. “XY or 2D?” therefore compares different properties.
The Hilbert-space dimension counts possible quantum states. Each site has four local occupations, so 36 sites have a full occupation space containing 436 basis vectors. This is a third meaning of dimension. “Full dimensional” is not a clear model name here; use “full 2D Fermi–Hubbard model”.
2. One lattice site: four states and two qubits
We retain one spatial orbital per site. An electron in it can have spin up (↑) or spin down (↓). A mode is a site-and-spin combination. The Pauli exclusion principle allows at most one electron per mode. Two electrons with opposite spins can therefore occupy the same site.
|0\rangle_i\leftrightarrow|00\rangle,\quad |\uparrow\rangle_i\leftrightarrow|10\rangle,\quad |\downarrow\rangle_i\leftrightarrow|01\rangle,\quad |\uparrow\downarrow\rangle_i\leftrightarrow|11\rangle.
\]
In this local ket notation the up occupation is written first and the down occupation second. Qubit value 1 means “this mode is occupied”; 0 means “empty”. The qubits encode occupation. A qubit is not simply the spin arrow of one permanently assigned electron. Thus 36 sites require 72 modes and, in this direct encoding, 72 qubits even when only 32 electrons are present.
\dim\mathcal H=4^{36}=2^{72},\qquad \dim\mathcal H_{N_\uparrow=N_\downarrow=16}=\binom{36}{16}^{2}.
\]
The second expression counts only states with exactly sixteen electrons of each spin. The dynamics stays in that sector. A large state space alone does not establish quantum advantage: structure, entanglement, time and the requested observables also affect classical cost.
3. From occupation to operators
A quantum state is a superposition of occupation patterns. The complex coefficient an is an amplitude; its squared magnitude is the probability of measuring that pattern. An operator acts on a state. The Hamiltonian H is the energy operator and also determines how the state changes in time.
|\psi\rangle=\sum_{\mathbf n}a_{\mathbf n}|\mathbf n\rangle,\qquad \sum_{\mathbf n}|a_{\mathbf n}|^2=1.
\]
The operator c†iσ adds an electron with spin σ at site i; ciσ removes one. Adding to an already occupied mode gives zero. The number operator niσ counts that mode’s occupation:
n_{i\sigma}=c^\dagger_{i\sigma}c_{i\sigma},\qquad n_{i\sigma}|\mathbf n\rangle=n_{i\sigma}^{(\mathbf n)}|\mathbf n\rangle,\qquad n_{i\sigma}^{(\mathbf n)}\in\{0,1\}.
\]
Fermionic operators obey the anticommutation rules below. Here {A,B} = AB + BA, and δ equals 1 when the indices match and 0 otherwise. Exchanging two creation operators for different modes introduces a minus sign. That sign changes interference between amplitudes and must therefore be preserved in the simulation.
\{c_{i\sigma},c^\dagger_{j\sigma’}\}=\delta_{ij}\delta_{\sigma\sigma’},\qquad \{c_{i\sigma},c_{j\sigma’}\}=0.
\]
4. Hopping: what exactly does t do?
For two connected sites i and j, c†iσcjσ describes motion from j to i while preserving spin. Operators act from right to left: remove at j, then add at i. The reverse motion is the Hermitian-conjugate term. Together they make the energy operator Hermitian, giving real energy eigenvalues.
H_{ij,\sigma}=-t\left(c^\dagger_{i\sigma}c_{j\sigma}+c^\dagger_{j\sigma}c_{i\sigma}\right).
\]
t is an energy, not a time or a count of hops. Consider one electron on two sites. In the basis |L⟩, |R⟩, with the electron on the left or right, the Hamiltonian is a 2 × 2 matrix. Its energy eigenvalues are −t and +t. Starting on the left produces a coherent superposition:
H_{\rm one} = \begin{pmatrix}0&-t\\-t&0\end{pmatrix},\qquad |\psi(s)\rangle=\cos(ts/\hbar)|L\rangle+\mathrm{i}\sin(ts/\hbar)|R\rangle.
\]
P_R(s)=\sin^2(ts/\hbar).
\]
Here s denotes physical time, leaving t exclusively for the hopping energy. The ratio t/ℏ sets a frequency scale. The probability oscillates quantum mechanically; t is not a classical hopping rate.
5. Local repulsion U and double occupation
The interaction at one site is U ni↑ni↓. The occupation product is 1 only when both spins are present. A doubly occupied site is also called a doublon. In the order empty, up, down, doubly occupied, the interaction matrix is:
H_{U,i}=U n_{i\uparrow}n_{i\downarrow}=\operatorname{diag}(0,0,0,U).
\]
U > 0 means repulsion: double occupation adds energy. For an occupation pattern with d doubly occupied sites, the interaction contribution is Ud. This is only the interaction energy; the total energy also contains hopping. In a superposition the mean doublon count need not be an integer.
\langle H_U\rangle=U\sum_i\langle n_{i\uparrow}n_{i\downarrow}\rangle.
\]
A large positive U makes double occupation energetically costly, but does not forbid it. Excluding doubly occupied states would require an additional projection or an appropriate infinite-U limit. At U/t = 8 they remain in the state space and can appear during the dynamics.
6. The full Hamiltonian, term by term
\begin{aligned}H&=H_t+H_{t’}+H_U,\\ H_t&=-t\sum_{\langle i,j\rangle,\sigma}\left(c^\dagger_{i\sigma}c_{j\sigma}+c^\dagger_{j\sigma}c_{i\sigma}\right),\\ H_{t’}&=-t’\sum_{\langle\!\langle i,j\rangle\!\rangle,\sigma}\left(c^\dagger_{i\sigma}c_{j\sigma}+c^\dagger_{j\sigma}c_{i\sigma}\right),\\ H_U&=U\sum_i n_{i\uparrow}n_{i\downarrow}.\end{aligned}
\]
The first sum covers horizontal and vertical neighbours; the second covers diagonal neighbours. Each unordered pair is counted once, with both spins included for every pair. An open 6 × 6 lattice has 60 horizontal/vertical bonds and 50 diagonal bonds, plus 36 local U terms. Thus the number 8 is a coefficient, not a count of these terms.
We choose t > 0. Since t′ = −0.25t, the coefficient −t′ in the Hamiltonian is +0.25t. Its relative sign matters for interference. This target uses μ = 0. The code also supports a uniform −μN term; within a fixed-N sector it contributes only a global phase to time evolution.
Literature for this Hamiltonian: Hubbard’s original paper [1] and the review by Arovas et al. [2].
7. U = 8 means U/t = 8: writing out the units
t, t′ and U all have units of energy. A statement such as “U = 8” is therefore incomplete until an energy unit is specified. Here we choose the positive hopping energy t as the unit and divide the entire Hamiltonian by t. Call the resulting dimensionless operator h:
h=\frac{H}{t},\qquad u=\frac{U}{t}=8,\qquad r=\frac{t’}{t}=-0.25.
\]
\begin{aligned}h={}&-\sum_{\langle i,j\rangle,\sigma}(c^\dagger_{i\sigma}c_{j\sigma}+c^\dagger_{j\sigma}c_{i\sigma})\\ &+0.25\sum_{\langle\!\langle i,j\rangle\!\rangle,\sigma}(c^\dagger_{i\sigma}c_{j\sigma}+c^\dagger_{j\sigma}c_{i\sigma})\\ &+8\sum_i n_{i\uparrow}n_{i\downarrow}.\end{aligned}
\]
“t = 1, U = 8” is a choice of units: one energy unit is t and the local double-occupation energy is eight such units. If a hypothetical physical realisation had t = 0.25 eV, then U would be 2 eV. At t = 0.10 eV the same ratio would give U = 0.80 eV. These examples do not identify our model with a real material.
For time, start with the Schrödinger equation. Let s be physical time, ℏ the reduced Planck constant and τ = ts/ℏ a dimensionless time. Substitution gives:
\mathrm{i}\hbar\frac{\partial|\psi\rangle}{\partial s}=H|\psi\rangle,\qquad \tau=\frac{ts}{\hbar}\quad\Longrightarrow\quad \mathrm{i}\frac{\partial|\psi\rangle}{\partial\tau}=h|\psi\rangle.
\]
|\psi(\tau)\rangle=e^{-\mathrm{i}h\tau}|\psi(0)\rangle,\qquad T=0.4\ \Longleftrightarrow\ s=0.4\frac{\hbar}{t}.
\]
The project notation T = 0.4 means the final value τ = 0.4. With the hypothetical t = 0.25 eV this corresponds to about 1.05 femtoseconds of model time. The quantum computer implements gates and repeated measurements to represent that evolution. Execution time, queueing time and analysis time are separate quantities. Without a physical value of t, no time in seconds is specified.
8. Why XY terms appear in a fermion model
The Pauli operators X, Y and Z act on the two-dimensional space of one occupation qubit. I is the identity. In the basis |0⟩, |1⟩ their matrices are:
X=\begin{pmatrix}0&1\\1&0\end{pmatrix},\qquad Y=\begin{pmatrix}0&-\mathrm{i}\\\mathrm{i}&0\end{pmatrix},\qquad Z=\begin{pmatrix}1&0\\0&-1\end{pmatrix}.
\]
Thus n = (I − Z)/2 is exactly 0 on |0⟩ and 1 on |1⟩. Also, (X + iY)/2 = |0⟩⟨1| removes an occupation. To reproduce the fermionic algebra as well, Jordan–Wigner adds Z operators over all preceding modes in a chosen ordering:
c_p=\left(\prod_{k\lt p}Z_k\right)\frac{X_p+\mathrm{i}Y_p}{2},\qquad c_p^\dagger=\left(\prod_{k\lt p}Z_k\right)\frac{X_p-\mathrm{i}Y_p}{2}.
\]
For p < q this gives the hopping identity below. Operators on different qubits act as tensor products; for example XpXq means X on p and X on q. The Z string between the endpoints is the parity string:
c_p^\dagger c_q+c_q^\dagger c_p=\frac12\left(X_p Z_{p+1}\cdots Z_{q-1}X_q+Y_p Z_{p+1}\cdots Z_{q-1}Y_q\right).
\]
\prod_{p\lt k\lt q}Z_k\,|\mathbf n\rangle=(-1)^{\sum_{p\lt k\lt q}n_k}|\mathbf n\rangle.
\]
An even number of occupied intervening modes gives +1; an odd number gives −1. For consecutive modes the string is empty and equals I. In that case fermionic hopping reduces exactly to an XY term:
H_{p,p+1}=-\frac{t}{2}(X_pX_{p+1}+Y_pY_{p+1}).
\]
This explains why a correct Hubbard simulation can contain XY gates. The reverse inference does not hold: a collection of bare XY couplings does not automatically contain the required fermionic signs, U terms and diagonal hopping. Whether a string is empty depends on the ordering of modes, not just lattice distance. Even a horizontal nearest-neighbour bond can have an intervening mode in a spin-interleaved ordering.
A step-by-step derivation with executable examples is available in the official OpenFermion documentation [3].
9. An example where bare XY gives the wrong sign
Take three ordered modes 0, 1 and 2 and move an electron between 0 and 2. Below their occupations are written from left to right as n₀n₁n₂. With the middle mode empty, hopping gives amplitude −t. With it occupied, fermionic parity changes this to +t:
\langle100|H_{02}|001\rangle=-t,\qquad \langle110|H_{02}|011\rangle=+t.
\]
The bare operator −t(X₀X₂ + Y₀Y₂)/2 gives −t in both cases. On an isolated bond this difference can remain hidden in an occupation probability; for interfering paths and superpositions it is generally physically relevant. Omitting the parity string is therefore not a general optimisation of the same 2D fermion model.
10. How U and fermionic swaps act on qubits
For the up and down modes of one site, the interaction follows directly from n = (I − Z)/2. It can be implemented as a controlled phase: only |11⟩ receives a phase. In dimensionless time Δτ that phase is −uΔτ:
U n_{i\uparrow}n_{i\downarrow}=\frac{U}{4}(I-Z_{i\uparrow}-Z_{i\downarrow}+Z_{i\uparrow}Z_{i\downarrow}).
\]
e^{-\mathrm{i}u\Delta\tau\,n_{i\uparrow}n_{i\downarrow}}=\operatorname{diag}(1,1,1,e^{-\mathrm{i}u\Delta\tau}).
\]
At u = 8 and Δτ = 0.2 the angle is −1.6 radians. If this local interaction acted alone, occupation probabilities would remain unchanged: it changes only relative phase. Combined with hopping, that phase affects later interference and therefore occupations.
An alternative correct implementation brings modes next to one another using fermionic swaps. An ordinary SWAP exchanges two qubits; an fSWAP additionally gives |11⟩ a minus sign. Fermionic parity can thus be incorporated into routing instead of repeatedly using an explicit long Z string:
\mathrm{fSWAP}=\mathrm{SWAP}\,\mathrm{CZ},\qquad |00\rangle\mapsto|00\rangle,\quad |01\rangle\leftrightarrow|10\rangle,\quad |11\rangle\mapsto-|11\rangle.
\]
The project code contains an explicit Jordan–Wigner route and fSWAP routes. Both can represent the same full Hamiltonian. Their time-step and gate ordering can nevertheless differ, producing different approximation errors. Model equivalence does not automatically mean equal circuit error.
For fermionic swap networks, see Kivlichan et al. [4].
11. Worked example: two sites and two electrons at U/t = 8
To see the interplay of t and U exactly, consider a small auxiliary model: two sites, two electrons and one hopping bond. There is no diagonal t′ bond here. This is a teaching example; its result is not a prediction for our 6 × 6 lattice or its initial state.
Choose the mode order 1↑, 1↓, 2↑, 2↓. Define the singlet state |S⟩, with one electron at each site, and the symmetric doublon state |D₊⟩. The creation operators below also specify the sign convention:
|S\rangle=\frac{c^\dagger_{1\uparrow}c^\dagger_{2\downarrow}-c^\dagger_{1\downarrow}c^\dagger_{2\uparrow}}{\sqrt2}|\mathrm{vac}\rangle,\qquad |D_+\rangle=\frac{c^\dagger_{1\uparrow}c^\dagger_{1\downarrow}+c^\dagger_{2\uparrow}c^\dagger_{2\downarrow}}{\sqrt2}|\mathrm{vac}\rangle.
\]
The interaction energy of |S⟩ is zero; that of |D₊⟩ is U. Hopping mixes the states with matrix element −2t. In this two-dimensional subspace:
H_{S,D_+}=\begin{pmatrix}0&-2t\\-2t&U\end{pmatrix},\qquad \det(H-EI)=E^2-UE-4t^2=0.
\]
E_\pm=\frac{U\pm\sqrt{U^2+16t^2}}{2},\qquad \frac{E_-}{t}\bigg|_{U/t=8}=\frac{8-\sqrt{80}}{2}\simeq-0.4721.
\]
The lowest energy is negative because of hopping-induced mixing even though the local repulsion is positive. The probability of one doublon in this ground state is the Hellmann–Feynman derivative ∂E₋/∂U, since ∂H/∂U is the doublon number operator:
\langle D_{\rm tot}\rangle=\frac{\partial E_-}{\partial U}=\frac12\left(1-\frac{U}{\sqrt{U^2+16t^2}}\right)\simeq0.0528\quad(U/t=8).
\]
There is therefore about a 5.28% probability of double occupation somewhere in this dimer ground state. The mean per site is half that value. U = 8t strongly suppresses double occupation here but does not eliminate it. For U ≫ t, E₋ ≈ −4t²/U and the familiar superexchange scale J ≈ 4t²/U emerges. Such a spin description is a low-energy approximation; our full Hubbard dynamics also retains charge motion and doublon states.
For the strong-coupling limit and its relationship to spin models, see Arovas et al., §5.2 [2].
12. From the Hamiltonian to a finite quantum circuit
Exact evolution is exp(−ihT). Hopping and interaction terms generally do not commute: AB ≠ BA. Executing them separately in sequence is therefore an approximation. A symmetric product formula splits a time step into a forward and reverse order. For h = a + b the basic scheme is:
e^{-\mathrm{i}(a+b)\Delta\tau}=e^{-\mathrm{i}a\Delta\tau/2}e^{-\mathrm{i}b\Delta\tau}e^{-\mathrm{i}a\Delta\tau/2}+O(\Delta\tau^3),\qquad T=m\Delta\tau.
\]
The explicit reference route applies each hopping term for half a step, then the U terms for a full step, then the same hopping list in reverse order for half a step. The handbook specifies T = 0.4 with two steps of 0.2 for this target. T determines model time; the number of steps determines how that evolution is approximated. Other runs may use a different step subdivision.
At fixed T the usual asymptotic error order of this second-order formula is O(Δτ²). Its coefficient depends on commutators and model parameters, among other things. This provides no numerical error guarantee without further checks. Hardware errors add another contribution. A full Hamiltonian in the circuit means the correct model terms are included; accuracy separately requires time-step refinement and execution checks.
The dependence of product-formula errors on commutators is developed by Childs et al. [5].
13. What do we predict for the 6 × 6 lattice?
The initial state is an occupation pattern with alternating spins and four electrons removed relative to one electron per site. This leaves 32 electrons: sixteen up and sixteen down. “Four holes” means a deficit relative to 36 electrons, not relative to the maximum occupation of 72. The net hole doping is 4/36 ≈ 11.1%. During evolution additional empty sites can form together with doublons while N remains fixed.
N_\sigma=\sum_i n_{i\sigma},\qquad [H,N_\uparrow]=[H,N_\downarrow]=0,\qquad N_\uparrow=N_\downarrow=16.
\]
Conservation of N = 32 is known in advance. It is a useful execution check, but does not tell us where the electrons are. The local predictions are charge Cᵢ, spin imbalance Sᵢ and double-occupation probability Dᵢ. Brackets ⟨O⟩ mean the expectation value ⟨ψ(T)|O|ψ(T)⟩:
C_i=\langle n_{i\uparrow}+n_{i\downarrow}\rangle,\qquad S_i=\langle n_{i\uparrow}-n_{i\downarrow}\rangle,\qquad D_i=\langle n_{i\uparrow}n_{i\downarrow}\rangle.
\]
Here Sᵢ is a dimensionless spin imbalance. The physical z component of spin is ℏSᵢ/2. The allowed ranges are 0 ≤ Cᵢ ≤ 2, −1 ≤ Sᵢ ≤ 1 and 0 ≤ Dᵢ ≤ 1. Double occupation is a correlation: Dᵢ generally does not equal the product of the two mean occupations.
\begin{aligned}C_i&=1-\frac{\langle Z_{i\uparrow}\rangle+\langle Z_{i\downarrow}\rangle}{2},\\ S_i&=\frac{\langle Z_{i\downarrow}\rangle-\langle Z_{i\uparrow}\rangle}{2},\\ D_i&=\frac{1-\langle Z_{i\uparrow}\rangle-\langle Z_{i\downarrow}\rangle+\langle Z_{i\uparrow}Z_{i\downarrow}\rangle}{4}.\end{aligned}
\]
These formulas connect the model directly to measured qubit occupations. The question is whether the circuit predicts these local time-dependent quantities accurately enough. The model is inspired by correlated-electron physics. This short-time simulation does not establish superconductivity or an accurate description of a real cuprate material.
A related 72-qubit experiment is the 2D Hubbard study by Alam et al. on Google Willow [6]. It provides experimental context; model parameters, preparation and validation must be compared for each experiment.
Sources and verification of the model choice
Project conventions were checked in theory/HANDBOOK.md, cuprate_2d_hubbard_pipeline.py, cuprate_2d_fermionic_dynamics.py and cuprate_2d_fermionic_mapping.py in fermi-hubbard-2d-nighthawk. The explicit route carries the model label full_2d_fermi_hubbard and records fermionic parity, t′ and U. Small-lattice checks are in validate_cuprate_2d_fermionic_dynamics.py. This identifies the model and verification procedure; it does not automatically certify every large hardware run.
Literature and further reading
-
J. Hubbard. Electron correlations in narrow energy bands · Free full text. Proceedings of the Royal Society A 276, 238–257 (1963).
Original paper on a model of mobile electrons with local interactions. Historical foundation for the Hamiltonian in section 6; our particular lattice size and parameters are project choices. -
D. P. Arovas, E. Berg, S. A. Kivelson & S. Raghu. The Hubbard Model · arXiv. Annual Review of Condensed Matter Physics 13, 239–274 (2022).
Review of the Hubbard Hamiltonian, two-dimensional lattices and controlled limits. The introduction and §5.2 on strong coupling are particularly useful alongside the discussion of U/t and the superexchange scale in section 11. -
Google Quantum AI / OpenFermion. The Jordan–Wigner and Bravyi–Kitaev Transforms.
Official documentation covering anticommutation, occupation encoding and the Jordan–Wigner parity string. Includes executable examples relevant to the operators and XY relationship in sections 3, 8 and 9. -
I. D. Kivlichan et al.. Quantum Simulation of Electronic Structure with Linear Depth and Connectivity · arXiv. Physical Review Letters 120, 110501 (2018).
Research paper on fermionic swap networks: modes are brought next to one another while preserving fermionic signs. Supports the routing principle in section 10; our specific gate counts must come from our own circuits. -
A. M. Childs, Y. Su, M. C. Tran, N. Wiebe & S. Zhu. Theory of Trotter Error with Commutator Scaling · arXiv. Physical Review X 11, 011020 (2021).
Research paper on product-formula error bounds and their dependence on commutators. Further reading for the symmetric time steps in section 12. A numerical bound for our lattice requires applying the theory to our Hamiltonian and step size. -
F. Alam et al.. Programmable digital quantum simulation of 2D Fermi-Hubbard dynamics using 72 superconducting qubits. arXiv:2510.26845v2 (2025).
Experimental preprint on 2D Hubbard dynamics on Google Willow, reaching 6 × 6 sites with 72 qubits. Compare equation (1), the second-order Trotter approach and classical checks. Its Hamiltonian and initial conditions are not identical to our t–t′–U target; this paper does not validate our IBM results.
The literature supports the model and methods. Our specific settings U/t = 8, t′/t = −0.25 and T = 0.4, and the circuits actually executed, are recorded in the project sources listed above. The two-site examples are teaching derivations with their own assumptions.