Simulation of Lattice Models

Quantum simulation of a lattice model prepares, evolves, or measures a finite spin, fermion, boson, or constrained model on a quantum device. The graph, boundaries, filling, state, and observable specify the physical problem. Encoding and hardware placement then determine whether its local interactions remain inexpensive to implement.

This page develops that connection through a four-site Ising quench and contrasts spin, fermion, and boson registers. A finite-lattice result can test a mechanism; predicting a real substance additionally needs the model checks in Simulation of Quantum Materials.

Required background. What Is Quantum Simulation? defines the requested observable; Lattice Models Overview introduces local degrees of freedom; Hamiltonian Simulation and Digital Quantum Simulation develop evolution approximations and circuits.

Model, Encoding, and Hardware Locality

Locality is not merely a property of the formula written by the theorist. It must survive four representations:

1. model locality: each term acts on a bounded region of the lattice;
2. encoding locality: that term maps to an operator of manageable support on the register;
3. hardware locality: the required register interactions fit the device graph with limited routing;
4. observable locality: the requested output can be estimated without reconstructing an exponentially large state.

A nearest-neighbor fermion hopping term can become a long Pauli string after a poor mode ordering. A local spin interaction can require many SWAPs on a mismatched processor. Conversely, a globally defined quantity such as a structure factor may still be efficiently assembled from repeated local measurements. The correct resource estimate follows the entire chain, not the first Hamiltonian alone.

Finite Graphs, Boundaries, and Observables

Use $\hbar=1$ and dimensionless spin operators throughout, so the couplings have energy units and times have inverse-energy units. The finite model specifies a graph $G_L$, local spaces $\mathcal H_x$, Hamiltonian $H_L$, constraints and sector, initial state, and requested observables.

Geometry and support

Write the Hamiltonian as a sum over supports,

$$H_L = \sum_{X\subseteq G_L} h_X,$$

and record the sites on which each $h_X$ acts. For a finite-range model, one usually has

$$\operatorname{diam}(X)\le R, \qquad \lVert h_X\rVert\le J_0,$$

with interaction range $R$ and local scale $J_0$ independent of system size, and a bounded number of terms incident on each site. This support list, rather than a model nickname, determines parallel layers, commutators, routing, and oracle cost.

Boundaries are part of the Hamiltonian

Open, periodic, twisted, cylindrical, and irregular boundaries define different finite operators. For a nearest-neighbor ring,

$$H_{\mathrm{PBC}} = \sum_{j=1}^{L-1}h_{j,j+1} +h_{L,1},$$

whereas the open chain omits $h_{L,1}$. The wrapping term is geometrically local on the target ring even if it is nonlocal on a linear processor. This is a routing fact, not permission to omit the term. Boundary conventions and finite-size consequences are treated in Boundary Conditions on Lattices.

Sector and initial condition

State all exact charges used to restrict the calculation. Typical examples are total magnetization, particle number, fermion parity, momentum, and point group quantum numbers. Also distinguish:

$$\rho_0 = \lvert\psi_0\rangle\langle\psi_0\rvert, \qquad \rho_\beta = \frac{e^{-\beta H_L}}{Z_L},$$

where $\beta=1/(k_BT)$ and the thermal trace defining $Z_L$ is over the chosen sector. For a grand-canonical ensemble, use $H_L-\mu\hat N$ in the exponent and trace over particle sectors, counting the chemical-potential term only once. A pure-state quench and a thermal calculation are different tasks even when they use the same Hamiltonian.

Observable and normalization

A local expectation, extensive sum, and intensive density carry different precision requirements. If

$$O_{\mathrm{tot}}=\sum_{x=1}^{L}O_x, \qquad o_L=\frac{O_{\mathrm{tot}}}{L},$$

then fixed additive error in $O_{\mathrm{tot}}$ becomes increasingly strict per site, while fixed error in $o_L$ defines a stable bulk target. State which quantity is reported and whether the intended conclusion concerns this finite instance or a thermodynamic extrapolation.

Representative Model Families

The following Hamiltonians illustrate how local terms translate into register operations. Each edge is counted once; the hopping amplitudes below are real. For complex hopping, replace $t_{ij}(c_i^\dagger c_j+c_j^\dagger c_i)$ by $t_{ij}c_i^\dagger c_j+t_{ij}^*c_j^\dagger c_i$ (and similarly for bosons).

Family	Representative terms	Natural register issue	Typical outputs
Ising	$Z_iZ_j$, $X_i$, longitudinal fields	direct qubit map; edge scheduling	magnetization, domains, gaps, quenches
Heisenberg or XXZ	$X_iX_j+Y_iY_j+\Delta Z_iZ_j$	exchange compilation; spin symmetry	correlators, transport, structure factors
Fermi–Hubbard	$c_{i\sigma}^\dagger c_{j\sigma}+\mathrm{h.c.}$, $n_{i\uparrow}n_{i\downarrow}$	anticommutation, mode order, routing	densities, spin correlations, Green functions
Bose–Hubbard	$b_i^\dagger b_j+\mathrm{h.c.}$, $n_i(n_i-1)$	local occupation cutoff or native modes	number statistics, coherence, expansion

For example, a transverse-field Ising family can be written

$$H_{\mathrm I} = -\sum_{(i,j)\in E}J_{ij}Z_iZ_j -\sum_i h_iX_i.$$

The isotropic spin-$1/2$ Heisenberg family is

$$H_{\mathrm H} = \frac14 \sum_{(i,j)\in E}J_{ij} \left( X_iX_j+Y_iY_j+Z_iZ_j \right).$$

A spinful Hubbard instance is

$$H_{\mathrm F} ={} -\sum_{(i,j),\sigma} t_{ij} \left( c_{i\sigma}^\dagger c_{j\sigma} +c_{j\sigma}^\dagger c_{i\sigma} \right) + U\sum_i n_{i\uparrow}n_{i\downarrow} -\sum_{i,\sigma}\mu_{i\sigma}n_{i\sigma}.$$

A truncated Bose–Hubbard instance begins from

$$H_{\mathrm B} = -\sum_{(i,j)}t_{ij} \left( b_i^\dagger b_j+b_j^\dagger b_i \right) +\frac{U}{2}\sum_i n_i(n_i-1) -\sum_i\mu_i n_i.$$

The letters $L$, $E$, $J$, $t$, $U$, and $\mu$ are not enough to reproduce a calculation. One must still provide graph orientation, coefficient units, operator normalization, boundaries, local cutoffs, and sector data.

Encode the Local Degrees of Freedom

Spin-1/2 sites

One physical or logical qubit per site is the simplest exact mapping,

$$\lvert\uparrow\rangle_i\leftrightarrow\lvert0\rangle_i, \qquad \lvert\downarrow\rangle_i\leftrightarrow\lvert1\rangle_i, \qquad S_i^a=\frac12\sigma_i^a.$$

This agrees with $Z|0\rangle=|0\rangle$. Exchanging the two labels requires transforming the spin operators as well; under conjugation by $X$, $S^z$ and $S^y$ change sign while $S^x$ does not. A direct map preserves the support of every Pauli interaction. It does not guarantee low cost if the device coupling graph differs from the target graph.

Higher spins and qudits

A local dimension $q$ can use a native qudit, a binary encoding with $\lceil\log_2q\rceil$ qubits, or a unary encoding with $q$ basis states spread over $q$ qubits and a one-hot constraint. Binary encoding minimizes qubit count but can make local operators denser. Unary encoding uses more qubits but often turns transitions between adjacent levels into low-weight operations.

If $2^m>q$ for a binary register, the unused states form a leakage space. With $P$ the projector onto valid states, a simulation should report at least

$$p_{\mathrm{valid}} = \operatorname{Tr}(P\rho), \qquad p_{\mathrm{leak}}=1-p_{\mathrm{valid}}.$$

Rejecting invalid states changes the sampling cost. The accepted state is $P\rho P/p_{\mathrm{valid}}$ when $p_{\mathrm{valid}}>0$; it describes conditional outcomes, not the original unconditioned experiment.

Fermionic modes

Fermionic occupation has finite local dimension, but anticommutation requires parity signs. For an ordered set of modes and the convention that $\lvert1\rangle_j$ is occupied, Jordan–Wigner gives

$$c_j = \left( \prod_{k<j}Z_k \right) \frac{X_j+iY_j}{2}.$$

Consequently, for $p<q$, a hopping operator contains the intervening 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).$$

The derivation, occupation convention, and boundary signs are explained in Jordan–Wigner Transformation. For simulation, the key consequence is that mode order is a compilation choice. Serpentine order can keep many bonds short on a rectangular lattice, but no one-dimensional ordering keeps every two-dimensional nearest-neighbor bond adjacent. Auxiliary-qubit and locality-preserving encodings trade extra registers and constraints for shorter operator support.

Fermionic SWAP networks provide another strategy: move modes through one another while applying interactions when the relevant pair becomes adjacent. The logical layout then changes during the circuit, so the compiler must track which qubit currently represents each mode. An ordinary SWAP and a fermionic SWAP differ by the sign on the doubly occupied state.

Bosonic modes

A qubit register cannot exactly hold the infinite basis of an ideal bosonic mode. A local cutoff retains

$$\mathcal H_i^{(n_{\max})} = \operatorname{span} \left\{ \lvert0\rangle, \ldots, \lvert n_{\max}\rangle \right\},$$

with dimension $q=n_{\max}+1$. The truncated creation operator is

$$b_i^\dagger = \sum_{n=0}^{n_{\max}-1} \sqrt{n+1} \lvert n+1\rangle\langle n\rvert.$$

This is a new finite model, not an exact representation of the unbounded one. Convergence requires increasing $n_{\max}$ and monitoring upper-level weight, for example

$$p_{\mathrm{edge},i} = \langle n_{\max}\rvert\rho_i\lvert n_{\max}\rangle.$$

Small upper-level weight also need not bound virtual-excitation energy shifts; compare the observable at larger cutoffs. A small mean occupation alone is insufficient when rare high-occupation events control the observable. Native oscillator platforms can avoid a qubit cutoff at the register level, but finite anharmonicity, detector range, and leakage still impose an operational truncation.

Project a polynomial operator before replacing it by powers of truncated matrices. For $x=(b+b^\dagger)/\sqrt2$ and the occupation projector $P$, $Px^2P-(PxP)^2=(n_{\max}+1)|n_{\max}\rangle\langle n_{\max}|/2$. The missing term comes from a virtual step outside the retained space and back.

Symmetries and Constraints

An encoding is exact only on its declared physical space. Let $Q_a$ denote conserved charges and $G_x$ local constraints. The ideal target obeys

$$[H_L,Q_a]=0, \qquad G_x\lvert\psi_{\mathrm{phys}}\rangle =g_x\lvert\psi_{\mathrm{phys}}\rangle, \qquad (I-P_{\mathcal C})H_LP_{\mathcal C}=0.$$

The last condition makes the physical subspace invariant under evolution; merely identifying its basis states does not guarantee this. Simultaneous Pauli tapering additionally requires independent, mutually commuting generators. Symmetry can reduce qubit count, constrain an ansatz, enable tapering, and provide diagnostics. These are distinct uses. A circuit that commutes with particle number may still leave the intended momentum sector. Postselecting on a measured charge can mitigate detected leakage, but changes the accepted-shot count and does not repair coherent errors within the sector.

For an implemented state $\rho$, report sector leakage through projectors,

$$p_{\mathrm{out}} = 1-\operatorname{Tr}(P_{\mathcal C}\rho),$$

and condition observables only when the conditioning rule and its uncertainty were specified in advance. Symmetry Sectors develops the corresponding block structure for classical many-body methods.

Turn the Interaction Graph into Circuit Layers

The conflict graph

Suppose two-site terms are associated with target edges $e\in E$. Construct a conflict graph whose vertices are the terms and whose edges join terms that cannot execute simultaneously. A proper coloring partitions the Hamiltonian,

$$H_L = \sum_{a=1}^{\chi}H_a, \qquad H_a=\sum_{e\in\mathcal M_a}h_e,$$

where $\mathcal M_a$ is a matching when the only conflict is sharing a site. All gates in one matching can then run in parallel on hardware that supports the corresponding couplings.

For an open nearest-neighbor chain, even and odd bonds suffice:

$$\begin{aligned} H_{\mathrm{odd}} &= h_{1,2}+h_{3,4}+\cdots, \\ H_{\mathrm{even}} &= h_{2,3}+h_{4,5}+\cdots. \end{aligned}$$

An odd periodic ring needs three matching colors, while an even ring needs two (for simple cycle graphs with at least three sites). A $D$-dimensional hypercubic nearest-neighbor lattice has a $2D$ direction-and-parity schedule when boundaries are open or each periodic side has even length. An odd periodic side needs a third matching in that direction; the two-color wraparound schedule fails. These counts describe logical interaction depth; control constraints, crosstalk, and routing can increase physical depth.

Terms need not be disjoint to commute. For example, all Ising bonds $Z_iZ_j$ commute even when they share a site. They may be grouped into one exponential algebraically, yet still require several physical layers because one qubit cannot participate in two simultaneous two-qubit gates. Algebraic grouping and hardware scheduling are related but distinct.

Layout and routing

Let $G_T$ be the target interaction graph and $G_H$ the hardware coupling graph. A placement $\pi:V(G_T)\to V(G_H)$ is useful when most target edges map to short paths. A simple routing burden is

$$C_{\mathrm{route}}(\pi) = \sum_{(i,j)\in E(G_T)} w_{ij} \left[d_H(\pi(i),\pi(j))-1\right],$$

where $d_H$ is hardware graph distance and $w_{ij}$ weights how often the interaction is used. This is only a proxy: parallel congestion, directionality, native gate family, calibration quality, and dynamical remapping also matter.

The wraparound edge of a periodic chain illustrates the issue. It is local in $G_T$. On a linear $G_H$, its endpoints are far apart, so implementing it may require a SWAP network, teleportation, or a different mode layout. Calling the term “nonlocal” without naming the graph conflates model geometry with device geometry.

Track modes through swaps

A swap network changes which physical qubit represents each logical mode. Apply the final permutation to the measurement operators as well as to the Hamiltonian gates. Otherwise a correctly executed circuit can return the density of the wrong site. Fermionic swaps also supply the minus sign on double occupation; replacing them with ordinary swaps changes coherent many-fermion states (Kivlichan and colleagues, 2018).

Choose the Simulation Paradigm

Digital execution

A digital simulator represents the target with gates. It is attractive when the model family, boundaries, term strengths, or observables must change frequently. The implementation may use product formulas, randomized formulas, block-encoding methods, qubitization, variational circuits, or phase estimation. These choices have different ancilla, precision, and fault-tolerance requirements.

For a local decomposition $H=\sum_a H_a$, first-order evolution is

$$U_1(\delta t) = \prod_{a=1}^{\chi} e^{-iH_a\delta t},$$

and $r=t/\delta t$ steps approximate $e^{-iHt}$. The leading error depends on commutators between overlapping local terms, not simply on the square of the number of all terms. Locality-aware product-formula bounds can therefore be far tighter than a triangle-inequality estimate that treats every pair as noncommuting.

For fixed spatial dimension, finite interaction range, bounded coordination, bounded local term norms, and fixed finite local dimension, rigorous algorithms can achieve gate count nearly linear in spacetime volume, up to subpolynomial or polylogarithmic factors under their stated assumptions. This is an asymptotic result about encoded local Hamiltonians. State preparation, observable estimation, error correction, and a nonlocal encoding can still dominate an application. Increasing a bosonic cutoff changes local dimension and operator norms; a fixed-cutoff complexity statement is not uniform in that limit (Childs–Su, 2019; Haah and colleagues, 2023).

Analog execution

An analog simulator realizes a laboratory generator $H_{\mathrm{lab}}$ whose effective dynamics approximate a target family. Optical lattices naturally realize Hubbard-like motion and interactions; neutral-atom and ion arrays naturally realize programmable spin models. Analog execution can preserve spatial locality without translating every term into gates, but it replaces gate-synthesis error with calibration, unwanted couplings, leakage, and effective-Hamiltonian error.

The target claim needs an explicit map

$$P H_{\mathrm{lab}}P = \lambda H_L+\beta I+\Delta H,$$

where $P$ selects the encoded subspace, $\lambda$ fixes the time scale, $\beta$ is dynamically irrelevant only within that subspace, and $\Delta H$ contains residual terms. Analog Quantum Simulation derives this correspondence. The projected equality alone does not control $(I-P)H_{\mathrm{lab}}P$: population can leave the subspace, and weak virtual excursions can shift its energies. Bound leakage and effective energy corrections on the simulated time interval, or use an exactly invariant encoding. A scalar shift is observable if a coherent reference branch experiences a different shift.

Hybrid execution

Hybrid simulation uses a quantum device for a state, overlap, response kernel, or local evolution while a classical loop selects parameters, embeds a cluster, extrapolates size, or reconstructs a spectrum. Examples include variational ground states, quantum-assisted imaginary-time updates, and cluster or embedding workflows. The correct resource is the cost of the whole loop, including failed optimizer iterations and repeated quantum calls.

No paradigm is uniformly best. A native analog array may dominate for broad short-time correlation maps, a shallow digital circuit may be best for a small programmable instance, and fault-tolerant qubitization may win at stringent energy precision. The observable, precision, and available controls determine the useful route.

Prepare the Requested State

Product states and quenches

Computational-basis spin configurations, occupation bitstrings, charge-density waves, domain walls, and Néel states are often inexpensive. They are ideal for quench dynamics because preparation error can be characterized locally and short-time behavior can be checked analytically.

A quench specifies both Hamiltonians,

$$\lvert\psi(t)\rangle = e^{-iH_{\mathrm f}t} \lvert\psi_0(H_{\mathrm i})\rangle.$$

It is not enough to state the final couplings. The initial preparation, quench profile, clock origin, and any ramp time belong to the protocol.

Adiabatic preparation

Choose a path

$$H(s) =(1-s)H_0+sH_1, \qquad s=t/T,$$

or a problem-specific nonlinear schedule. Success depends on the path, minimum relevant gap, matrix elements of $\partial_sH$, degeneracy structure, and endpoint overlap. A large final gap does not rule out an avoided crossing along the path. At a finite-size critical bottleneck, the required time can grow rapidly with $L$.

Symmetry can help or hurt. Remaining in the desired exact sector can avoid irrelevant crossings, but a path that preserves the wrong symmetry cannot reach a target state in another sector. Report the schedule and verify convergence with $T$ rather than labeling a single ramp adiabatic.

Variational preparation

For a parameterized state $\lvert\psi(\boldsymbol\theta)\rangle$, ground-state preparation often minimizes

$$E(\boldsymbol\theta) = \langle\psi(\boldsymbol\theta)\rvert H_L \lvert\psi(\boldsymbol\theta)\rangle.$$

Hamiltonian variational ansätze built from lattice term groups can preserve symmetries and use local gates,

$$U(\boldsymbol\theta) = \prod_{p=1}^{P} \prod_{a=1}^{\chi} e^{-i\theta_{p,a}H_a}.$$

Low energy is not by itself proof of high fidelity when the gap is small or unknown. Optimization bias, estimator noise, ansatz restriction, and hardware error must be separated. VQE develops the variational algorithm.

Eigenstate filtering and phase estimation

Phase estimation can resolve eigenenergies from an input

$$\lvert\psi_{\mathrm{in}}\rangle = \sum_n c_n\lvert E_n\rangle,$$

but the probability of obtaining $E_n$ is $\lvert c_n\rvert^2$. High-precision evolution does not create overlap with the desired state. Resource estimates must include preparation and repetitions from imperfect overlap. See Quantum Phase Estimation.

Thermal states

A Gibbs target

$$\rho_\beta = \frac{e^{-\beta H_L}}{\operatorname{Tr}e^{-\beta H_L}}$$

can require purification, imaginary-time methods, probabilistic filtering, open-system engineering, or sampling from a variational ensemble. Low temperature combines ground-state preparation difficulty with entropy and normalization estimation. A state whose energy density matches a thermal value need not reproduce the Gibbs ensemble, especially in integrable or many-body localized regimes.

Evolve Without Losing the Declared Physics

Product-formula structure

For two groups $A+B=H$, the symmetric second-order step is

$$S_2(\delta t) = e^{-iA\delta t/2} e^{-iB\delta t} e^{-iA\delta t/2}.$$

Over fixed total time, its leading global error scales as $O(t\delta t^2)$ times norms of nested commutators under standard bounded-operator assumptions. The constant and size dependence matter. A convergence study should vary $\delta t$, term ordering, and preferably formula order while holding the physical task fixed.

Product formulas can preserve a charge exactly when every compiled group commutes with it:

$$[H_a,Q]=0 \quad\text{for every }a \quad\Longrightarrow\quad [U_1(\delta t),Q]=0.$$

This is stronger than conservation only in the $\delta t\to0$ limit and is a good reason to choose symmetry-respecting groupings.

Local versus global accuracy

Operator-norm accuracy for the entire $L$-site unitary is a strong guarantee. Many scientific tasks ask only for a local observable over finite time. A Lieb–Robinson bound limits the influence of distant perturbations for short-range systems, so an observable-specific simulation can sometimes use a finite causal neighborhood or tolerate errors far outside it. Its initial reduced state must also be available or prepared with a controlled error. This does not justify discarding distant terms arbitrarily: the truncation radius and resulting error must be derived for the interaction decay, time, and observable.

Long-range interactions

Power-law interactions,

$$J_{ij}\sim\frac{J_0}{r_{ij}^{\gamma}},$$

change term count, scheduling, causal bounds, and analog-model discrepancy. Truncating at radius $R_c$ defines

$$\Delta H_{R_c} = \sum_{r_{ij}>R_c}J_{ij}O_iO_j.$$

Its operator norm can grow with system size even when each omitted coupling is small. Validate the final observable against several cutoffs or include the tail in the target. Never relabel an uncontrolled hardware tail as harmless because the desired textbook model is nearest-neighbor.

Measure the Scientific Output

Local observables and correlations

Single-site polarization and density are

$$m_i^a(t)=\langle\sigma_i^a(t)\rangle, \qquad n_i(t)=\langle\hat n_i(t)\rangle.$$

Connected equal-time correlations remove products of one-point means,

$$C_{ij}^{ab} = \langle O_i^aO_j^b\rangle -\langle O_i^a\rangle \langle O_j^b\rangle.$$

Each bitstring provides one joint sample for every same-basis $Z$ or number correlator; repeated preparations estimate their expectations, but noncommuting bases require additional settings. Covariance between estimators must be retained when many correlators are combined.

Structure factors

For positions $\mathbf r_j$, a static structure factor can be normalized as

$$S^{aa}(\mathbf q) = \frac1L \sum_{i,j} e^{i\mathbf q\cdot(\mathbf r_i-\mathbf r_j)} C_{ij}^{aa}.$$

State whether connected or raw correlations are used and whether $1/L$ or $1/L^2$ normalization is intended. Structure Factors explains the physical interpretation and normalization choices.

Dynamical response and spectra

A real-time correlator has the form

$$C_{AB}(t) = \operatorname{Tr} \left[ \rho\,A(t)B(0) \right].$$

Ancilla interferometry or controlled insertions, with suitable decompositions for nonunitary $A,B$, can recover its real and imaginary parts. Weak-probe linear response instead directly gives the retarded commutator; it does not in general determine the full unsymmetrized correlator without additional state information, such as equilibrium fluctuation–dissipation relations. The finite record is windowed before Fourier transformation,

$$\widetilde C_{AB}(\omega) = \int_{-T}^{T} dt\, w(t)e^{i\omega t}C_{AB}(t).$$

Thus frequency resolution, broadening, maximum reliable evolution time, sampling cadence, and window choice are part of the answer. A smooth spectrum can reflect the window rather than intrinsic lifetime. Fermionic Green functions additionally require parity-aware insertion of creation and annihilation operators. Spectral Functions defines the many-body spectra.

Finite systems and order parameters

On a finite system with an exact symmetry, an order parameter odd under that symmetry can have zero expectation in every symmetry eigenstate. For $M=\sum_i Z_i$ and a spin-flip symmetry $P$ satisfying $PMP^\dagger=-M$,

$$P\lvert\psi\rangle=p\lvert\psi\rangle \quad\Longrightarrow\quad \langle M\rangle=0.$$

One instead studies $\langle M^2\rangle$, long-distance correlations, susceptibility, a small symmetry-breaking field with a stated order of limits, or the low-lying finite-size tower. One cluster and one order parameter value do not establish a thermodynamic phase.

Worked Example: Four-Site Ising Quench

Consider an open four-site chain,

$$H = H_Z+H_X,$$

with

$$H_Z = -J \left( Z_1Z_2+Z_2Z_3+Z_3Z_4 \right), \qquad H_X=-h\sum_{j=1}^{4}X_j.$$

Take $\lvert\psi_0\rangle=\lvert0000\rangle$ and the convention $Z\lvert0\rangle=\lvert0\rangle$. The task is to estimate

$$m_z(t) = \frac14\sum_{j=1}^{4} \langle Z_j(t)\rangle$$

and the connected end-to-end correlation $C_{14}^{zz}(t)$.

Parallel interaction layers

All three $ZZ$ terms commute, but on two-qubit hardware they use two matchings,

$$H_A=-J(Z_1Z_2+Z_3Z_4), \qquad H_B=-JZ_2Z_3.$$

The exact $Z$-interaction exponential is

$$e^{-iH_Z\delta t} = e^{-iH_A\delta t} e^{-iH_B\delta t}.$$

The factors are algebraically order-independent, while the physical schedule has one parallel layer for bonds $(1,2)$ and $(3,4)$ and one layer for bond $(2,3)$.

Symmetric time step

Use

$$S_2(\delta t) = e^{-iH_X\delta t/2} e^{-iH_Z\delta t} e^{-iH_X\delta t/2}.$$

With

$$R_X(\theta)=e^{-i\theta X/2}, \qquad R_{ZZ}(\theta)=e^{-i\theta Z\otimes Z/2},$$

one step uses $R_X(-h\delta t)$ on every site in each half layer and $R_{ZZ}(-2J\delta t)$ on every bond. For $r$ repeated steps, adjacent $X$ half layers merge. The ideal native-gate schedule therefore has $3r$ two-qubit gates, $2r$ two-qubit depth, and $r+1$ global $X$ layers. A compiler that reports only $3r$ entangling gates but ignores depth has omitted a key resource.

Where the formula error comes from

Using $[X,Z]=-2iY$,

$$[H_X,H_Z] = -2iJh \sum_{j=1}^{3} \left( Y_jZ_{j+1}+Z_jY_{j+1} \right).$$

The commutator vanishes when $J=0$ or $h=0$, so the product formula is exact in both solvable limits. A worst-case triangle bound is

$$\lVert[H_X,H_Z]\rVert \le 12\lvert Jh\rvert,$$

but an observable-level error can be smaller. The defensible procedure is to compare several $\delta t$ values against exact diagonalization for this small instance before extrapolating the schedule to larger $L$.

Checks with analytic limits

If $h=0$, the input is an eigenstate and

$$m_z(t)=1, \qquad C_{14}^{zz}(t)=0.$$

If $J=0$, sites rotate independently,

$$m_z(t)=\cos(2ht), \qquad C_{14}^{zz}(t)=0.$$

At general $J/h$, monitor norm, reflection symmetry, the measured energy after the quench, and the $\delta t^2$ convergence expected from the symmetric formula in the asymptotic regime. If periodic boundaries are added, the $Z_4Z_1$ bond changes both the target Hamiltonian and the hardware schedule.

Error Budget at the Observable Level

To add contributions, first bound their effects on the same observable in the same units. For a sequence of intermediate calculations $o_j$, bounds $|o_{j+1}-o_j|\le\epsilon_j$ give $|o_{\mathrm{final}}-o_0|\le\sum_j\epsilon_j$ by the triangle inequality. Correlation does not invalidate this deterministic bound, but standard deviations cannot simply be substituted for its terms. If component bounds are probabilistic, their simultaneous coverage also needs a failure-probability allocation (for example, a union bound).

Contribution	What it means	Required check
finite size and boundary	difference between finite target and bulk claim	several sizes, shapes, or twists
local cutoff	omitted boson or qudit levels	increase cutoff; monitor boundary weight
preparation	input differs from declared state or ensemble	state observables, energy, symmetry, overlap bounds
algorithm	formula, polynomial, filtering, optimizer, or time-window error	convergence in the controlling parameter
compilation	synthesis, routing, mode-tracking, or pulse approximation	circuit identities and independent compilation
hardware	noise, drift, crosstalk, leakage, and SPAM	interleaved calibration and raw data
statistics	finite accepted samples and covariance	confidence intervals and stopping rule
inference	fit, Fourier transform, mitigation, or extrapolation	sensitivity and held-out tests

Finite-size difference is not an error when the finite system itself is the declared target. It becomes an error only when the result is interpreted as a bulk property.

If an ideal output $\rho$ and implemented output $\widetilde\rho$ obey

$$D(\rho,\widetilde\rho) = \frac12\lVert\rho-\widetilde\rho\rVert_1 \le\delta,$$

then any bounded observable satisfies

$$\left| \operatorname{Tr} \left[O(\rho-\widetilde\rho)\right] \right| \le 2\lVert O\rVert_\infty\delta.$$

This converts a state-level guarantee into an observable bound, but it can be too pessimistic for local observables and too expensive to establish at large size. Observable-specific convergence and validation are therefore central.

Sampling and Derived Quantities

For independent samples $X_k$ of a bounded estimator with variance $\sigma_X^2$, the sample mean has

$$\operatorname{Var}(\overline X) = \frac{\sigma_X^2}{N}.$$

Reaching standard error $\epsilon_{\mathrm{stat}}$ therefore costs roughly $N\sim\sigma_X^2/\epsilon_{\mathrm{stat}}^2$, before confidence corrections. Autocorrelation, drift, rejected shots, and mitigation can reduce the effective sample size.

For a derived quantity $f(\boldsymbol\mu)$ with covariance matrix $\Sigma$, linear propagation gives

$$\operatorname{Var}(f) \approx \nabla f^{\mathsf T} \Sigma \nabla f.$$

This matters for connected correlations, structure factors, gaps obtained by subtraction, response derivatives, and fitted critical parameters. Treating correlated estimates as independent can either inflate or understate the final uncertainty.

Finite-Size Reasoning

Correlation length and boundaries

Away from criticality, bulk finite-size corrections often become small when the linear size exceeds the correlation length and measurements are far from open edges. Near a critical point, $\xi$ grows and this criterion fails. A finite-size scaling ansatz may take the form

$$O(L,g) = L^{-x_O} F_O\!\left((g-g_c)L^{1/\nu}\right) +L^{-x_O-\omega}G_O+\cdots.$$

The exponents, correction terms, aspect ratio, boundary condition, and fitting window are hypotheses to test, not decorations added after data collection. Thermodynamic Limit explains the order of limits.

Compare matched quantities

Do not compare the total energy of one size with energy density at another. For open systems, distinguish center observables from spatial averages, which mix bulk and boundary regions. For fermions, compare the same filling and sector; for twisted boundaries, compare the same flux convention. Odd–even and commensurability effects can exceed the quoted hardware uncertainty.

Verification and Validation

Different checks probe different possible errors. Retain tests that remain feasible as the system grows beyond exact diagonalization:

1. Operator tests: verify signs, normalizations, mode order, boundaries, and matrix elements on one- and two-site bases.
2. Exact small instances: compare complete distributions, time traces, and correlators against exact diagonalization.
3. Solvable limits: set hopping, interaction, or field to zero; use commuting and free-particle limits.
4. Short-time moments: compare derivatives generated by nested commutators, such as $d\langle O\rangle/dt=i\langle[H,O]\rangle$.
5. Convergence: vary step size, circuit depth, cutoff, ramp time, time window, and mitigation strength.
6. Conservation and leakage: track charges, constraints, normalization, valid-subspace probability, and energy when it should be conserved.
7. Independent formulations: change mapping, mode order, compiler, platform, or classical method while preserving the physical task.
8. Held-out predictions: reserve observables, times, sizes, or parameter points that were not used for calibration or fitting.

For a variational state, the energy variance

$$\operatorname{Var}_\psi(H) = \langle H^2\rangle_\psi -\langle H\rangle_\psi^2$$

vanishes for any state supported within one exact eigenspace. A small value means that the measured energy distribution is narrow: for $a>0$,

$$\Pr\{|E-\langle H\rangle|\ge a\} \le\frac{\operatorname{Var}_\psi(H)}{a^2}.$$

It does not imply proximity to a single eigenstate when several levels lie in that narrow interval, nor identify the ground state. Spectral gaps or additional state-sensitive measurements are needed. Variance-based experiments therefore combine it with symmetry and independent physical checks (Kokail and colleagues, 2019).

Resource Accounting

The first register estimate is

Target	Baseline logical storage
$L$ spin-$1/2$ sites	$L$ qubits
$L$ local $q$-level sites, binary	$L\lceil\log_2q\rceil$ qubits
$M$ finite fermionic modes	commonly $M$ qubits before reductions
$L$ bosonic sites with cutoff $n_{\max}$	$L\lceil\log_2(n_{\max}+1)\rceil$ qubits in binary

This table omits ancillas, gauge or locality auxiliaries, syndrome qubits, magic-state factories, and workspace for oracles. It also says nothing about depth or repetitions.

A useful total-work model is

$$N_{\mathrm{exec}} = N_{L} N_{g} N_{t} N_{\mathrm{prep}} N_{\mathrm{basis}} N_{\mathrm{noise}} N_{\mathrm{shot}},$$

where the factors count sizes, couplings, times, preparation variants, measurement settings, mitigation or calibration scales, and shots. Actual experiments may not factorize this neatly, but the expression exposes omitted sweeps. End-to-end cost should also include rejected runs, optimizer calls, classical preprocessing, routing, decoding, and postprocessing.

For fault-tolerant estimates, report logical qubits, non-Clifford count and depth, logical cycle time, code distance assumptions, factory footprint, failure budget, and repetitions. Resource Estimation Tools develops the full accounting workflow.

Choosing among Common Routes

Task feature	Often suitable route	Main qualification
short-time local dynamics from a product state	native analog or shallow product formula	generator calibration and reachable time
programmable finite spin model	direct qubit encoding with colored interaction layers	hardware graph and crosstalk
fermion lattice on qubits	parity-aware map plus optimized ordering or swap network	strings, constraints, and mode tracking
soft-core bosons	native modes or converged finite cutoff	occupation-tail evidence
ground energy with useful overlap	variational preparation or phase estimation	optimizer bias or fault-tolerant depth
many local observables	shared bitstrings, grouped bases, or randomized measurements	covariance and postprocessing cost
high-resolution spectrum	long coherent evolution plus controlled reconstruction	time window and state overlap
bulk phase claim	several sizes, boundaries, and matched observables	extrapolation dominates one-device precision

Classical competition must be task-specific. One-dimensional gapped ground states, free limits, stabilizer circuits, low-entanglement short-time dynamics, sign-problem-free Monte Carlo regimes, and small clusters can be classically tractable even when the full model family is hard. A quantum experiment can be scientifically valuable without demonstrating computational advantage.

Exercises

1. Color open and periodic chains

Find the minimum number of disjoint-bond layers for an open chain with $L\ge2$ sites, an even periodic ring with $L\ge4$, and an odd periodic ring with $L\ge3$. Explain the odd-ring answer.

Solution

An open chain needs two layers when $L\ge3$: color alternating bonds odd and even. For $L=2$, one layer suffices. An even ring also needs two because its cycle graph is bipartite and alternating edge colors close consistently.

An odd ring needs three. Alternating two colors around an odd cycle returns to the first vertex with the same color on both incident edges, which is forbidden. A third color resolves the conflict. This is a logical matching count; a hardware graph without the wrapping edge can require additional routing depth.

2. Count the four-site Ising circuit

For the worked example, use $r$ symmetric steps and assume native $R_{ZZ}$ gates on a line. Count the two-qubit gates, two-qubit depth, and global $X$-rotation layers after merging adjacent half layers. How does a periodic $Z_4Z_1$ term change the answer on line and ring hardware?

Solution

Each step applies three $R_{ZZ}$ gates. Bonds $(1,2)$ and $(3,4)$ are parallel, followed by $(2,3)$, so the totals are

$$N_{2q}=3r, \qquad D_{2q}=2r.$$

The first and last steps retain an $X$ half layer, while every adjacent pair of half layers merges into a full layer. There are therefore $r+1$ global $X$-rotation layers.

On ring hardware, $(4,1)$ joins $(2,3)$ in the second matching, so there are $4r$ two-qubit gates but still depth $2r$. On line hardware, $(4,1)$ is not a native edge. The logical count is still four interactions per step, but SWAP or another nonlocal primitive adds gates and depth; the exact overhead depends on whether the layout is restored after each step.

3. Derive the Ising commutator

For an open $L$-site chain with

$$H_Z=-J\sum_{j=1}^{L-1}Z_jZ_{j+1}, \qquad H_X=-h\sum_{j=1}^{L}X_j,$$

derive $[H_X,H_Z]$ and a triangle-inequality norm bound. In which parameter limits is every product formula between $H_X$ and $H_Z$ exact?

Solution

Only $X_j$ operators on an endpoint of a bond fail to commute with that bond. Using $[X,Z]=-2iY$ gives

$$[H_X,H_Z] = -2iJh \sum_{j=1}^{L-1} \left( Y_jZ_{j+1}+Z_jY_{j+1} \right).$$

Every Pauli string has norm one, so

$$\lVert[H_X,H_Z]\rVert \le 4\lvert Jh\rvert(L-1).$$

If $J=0$ or $h=0$, one Hamiltonian group vanishes and the split evolution is exact at any step size. The bound also shows why a generic all-pairs estimate is unnecessarily pessimistic: only overlapping local terms contribute.

4. Compare fermionic mode orders

Place six spinless modes on a $2\times3$ rectangular lattice. Compare the Jordan–Wigner parity-string lengths for horizontal and vertical nearest- neighbor hoppings using row-major order

$$(1,1),(1,2),(1,3),(2,1),(2,2),(2,3).$$

Which target edges remain adjacent in the mode order?

Solution

Every horizontal edge connects consecutive mode indices, so there are no intervening $Z$ operators. The vertical edges connect index pairs $(1,4)$, $(2,5)$, and $(3,6)$. Each pair has two intervening modes and therefore two $Z$ operators in the parity string.

Thus four horizontal edges remain adjacent, while none of the three vertical edges does. A serpentine order changes which vertical edges are short but cannot make every edge of this two-dimensional graph adjacent in one linear ordering. The exercise concerns encoding support; hardware routing can add a second, independent cost.

5. Derive the truncated boson commutator

For

$$b^\dagger = \sum_{n=0}^{n_{\max}-1} \sqrt{n+1}\lvert n+1\rangle\langle n\rvert,$$

show that

$$[b,b^\dagger] = I-(n_{\max}+1) \lvert n_{\max}\rangle\langle n_{\max}\rvert.$$

What does its expectation value say about a cutoff test?

Solution

Direct multiplication gives

$$bb^\dagger = \sum_{n=0}^{n_{\max}-1} (n+1)\lvert n\rangle\langle n\rvert,$$

and

$$b^\dagger b = \sum_{n=1}^{n_{\max}} n\lvert n\rangle\langle n\rvert.$$

Subtracting leaves unit coefficient on every state below the cutoff and $-n_{\max}$ on the top state, which is the claimed operator. For a state with top-level probability $p_{\mathrm{edge}}$,

$$\langle[b,b^\dagger]\rangle = 1-(n_{\max}+1)p_{\mathrm{edge}}.$$

Small $p_{\mathrm{edge}}$ is a useful diagnostic, but convergence of the final observable under increasing $n_{\max}$ is still required.

6. Finite symmetry and magnetization

Let $P=\prod_i X_i$ be the global spin flip and $M=\sum_i Z_i$. Show that a parity eigenstate has $\langle M\rangle=0$ but can have nonzero $\langle M^2\rangle$. Evaluate both quantities in

$$\lvert\mathrm{cat}_+\rangle = \frac{ \lvert0\cdots0\rangle +\lvert1\cdots1\rangle }{\sqrt2}.$$

Solution

Because $PMP^\dagger=-M$ and $P\lvert\psi\rangle=p\lvert\psi\rangle$,

$$\langle M\rangle = \langle\psi\rvert P^\dagger(PMP^\dagger)P\lvert\psi\rangle = -\langle M\rangle,$$

so $\langle M\rangle=0$. The operator $M^2$ is even under $P$ and need not vanish. For the cat state, the two components have magnetizations $+L$ and $-L$, hence

$$\langle M\rangle=0, \qquad \langle M^2\rangle=L^2.$$

This is why a zero one-point order parameter on a finite symmetry eigenstate does not rule out strong ordered correlations.

7. Propagate correlation uncertainty

Estimate

$$C_{ij}=\mu_{ij}-\mu_i\mu_j$$

from jointly estimated means $\boldsymbol\mu=(\mu_{ij},\mu_i,\mu_j)$. Write the linearized variance in terms of their covariance matrix $\Sigma$.

Solution

The gradient is

$$\nabla C = \begin{pmatrix} 1\\ -\mu_j\\ -\mu_i \end{pmatrix}.$$

Therefore

$$\operatorname{Var}(C) \approx \begin{pmatrix} 1&-\mu_j&-\mu_i \end{pmatrix} \Sigma \begin{pmatrix} 1\\ -\mu_j\\ -\mu_i \end{pmatrix}.$$

When all three quantities come from the same bitstrings, the off-diagonal covariances are generally nonzero. Dropping them is not conservative in general because their contribution can have either sign.

There is also a finite-sample bias. With $N>1$ independent shared shots, $\widehat C=\overline{O_iO_j}-\overline O_i\,\overline O_j$ has expectation $(1-1/N)C$. Multiplying by $N/(N-1)$ gives the usual unbiased sample covariance; What Is Quantum Simulation? derives the shared-shot correction. Bias and variance are separate questions.

8. Small energy variance without a selected eigenstate

Let $H=\operatorname{diag}(0,\delta,1)$, where $0<\delta<1$, and take $|\psi\rangle=(|0\rangle+|1\rangle)/\sqrt2$. Find its energy, energy variance, and ground-state overlap. What happens as $\delta\to0^+$? Compare with the pure excited state $|1\rangle$.

Solution

For the superposition,

$$\langle H\rangle=\delta/2,\qquad \langle H^2\rangle=\delta^2/2,\qquad \operatorname{Var}(H)=\delta^2/4, \qquad |\langle0|\psi\rangle|^2=1/2.$$

Both energy error and variance tend to zero, while the ground-state overlap remains one half for every positive $\delta$. At exactly zero, the ground space becomes two-dimensional, so overlap with the ground-space projector is a different question. The state $|1\rangle$ has zero variance at every $\delta$ and is excited when $\delta>0$. Variance diagnoses energy spread, not state identity; a spectral gap sets the scale needed to turn an energy bound into a ground-state-overlap bound.

References

• I. Bloch, J. Dalibard, and S. Nascimbène, “Quantum simulations with ultracold quantum gases,” Nature Physics 8, 267–276 (2012), doi:10.1038/nphys2259.
• S. Bravyi and A. Kitaev, “Fermionic quantum computation,” Annals of Physics 298, 210–226 (2002), doi:10.1006/aphy.2002.6254.
• A. Browaeys and T. Lahaye, “Many-body physics with individually controlled Rydberg atoms,” Nature Physics 16, 132–142 (2020), doi:10.1038/s41567-019-0733-z.
• C. Cade, L. Mineh, A. Montanaro, and S. Stanisic, “Strategies for solving the Fermi–Hubbard model on near-term quantum computers,” Physical Review B 102, 235122 (2020), doi:10.1103/PhysRevB.102.235122.
• A. M. Childs and Y. Su, “Nearly optimal lattice simulation by product formulas,” Physical Review Letters 123, 050503 (2019), doi:10.1103/PhysRevLett.123.050503.
• A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, “Theory of Trotter error with commutator scaling,” Physical Review X 11, 011020 (2021), doi:10.1103/PhysRevX.11.011020.
• A. Elben et al., “Cross-platform verification of intermediate scale quantum devices,” Physical Review Letters 124, 010504 (2020), doi:10.1103/PhysRevLett.124.010504.
• I. M. Georgescu, S. Ashhab, and F. Nori, “Quantum simulation,” Reviews of Modern Physics 86, 153–185 (2014), doi:10.1103/RevModPhys.86.153.
• J. Haah, M. B. Hastings, R. Kothari, and G. H. Low, “Quantum algorithm for simulating real time evolution of lattice Hamiltonians,” SIAM Journal on Computing 52, FOCS18-250–FOCS18-284 (2023), doi:10.1137/18M1231511.
• I. D. Kivlichan, J. McClean, N. Wiebe, C. Gidney, A. Aspuru-Guzik, G. K.-L. Chan, and R. Babbush, “Quantum simulation of electronic structure with linear depth and connectivity,” Physical Review Letters 120, 110501 (2018), doi:10.1103/PhysRevLett.120.110501.
• C. Kokail et al., “Self-verifying variational quantum simulation of lattice models,” Nature 569, 355–360 (2019), doi:10.1038/s41586-019-1177-4.
• E. H. Lieb and D. W. Robinson, “The finite group velocity of quantum spin systems,” Communications in Mathematical Physics 28, 251–257 (1972), doi:10.1007/BF01645779.
• C. Monroe et al., “Programmable quantum simulations of spin systems with trapped ions,” Reviews of Modern Physics 93, 025001 (2021), doi:10.1103/RevModPhys.93.025001.
• J. D. Whitfield, V. Havlíček, and M. Troyer, “Local spin operators for fermion simulations,” Physical Review A 94, 030301(R) (2016), doi:10.1103/PhysRevA.94.030301.