VQE for H2 Excited States — Variational Quantum Deflation (VQD)

What You'll Learn:

• How to find excited states of molecules on a quantum computer, not just the ground state
• How the VQD algorithm uses overlap penalties to enforce orthogonality between solutions
• Why the β (penalty coefficient) must exceed the energy gap — and what happens when it doesn't
• How to compute state overlaps via statevector simulation (and how swap tests work on hardware)

Level: Advanced | Time: 30 minutes | Qubits: 2 | Framework: Qiskit

Prerequisites

• VQE for H2 Ground State — standard VQE algorithm, Pauli measurement
• Bell State — entanglement with CX gates
• Hadamard Test — expectation value estimation

The Idea

Standard VQE is like searching for the lowest point in a valley — the ground state. But chemistry doesn't stop at ground states. Photosynthesis, vision, fluorescence, and most spectroscopy depend on excited states — the higher-energy configurations a molecule can occupy.

VQD (Variational Quantum Deflation) extends VQE with a clever trick: after finding the ground state, it searches for the next state while penalizing overlap with the ground state. Think of it like a ball rolling downhill, but with invisible walls placed at every valley it has already found. The ball is forced to settle into the next valley — the first excited state.

For H2, the ground state is the bonding orbital (electrons shared between the two protons, a superposition of |01⟩ and |10⟩). The first excited state (a 2-fold degenerate pair involving |00⟩ and |11⟩ components) sits ~0.85 Ha above. The energy gap determines the molecule's absorption spectrum — the wavelength of light it absorbs.

How It Works

Step 1: Find the Ground State (Standard VQE)

This is identical to the H2 VQE tutorial. The optimizer minimizes:

L₀(θ) = ⟨ψ(θ)|H|ψ(θ)⟩

The variational principle guarantees convergence to the ground state energy E₀ ≈ -1.137 Ha.

Step 2: Find the First Excited State (VQD)

Now we run a second optimization with a modified cost function:

CODECopy
L₁(θ) = ⟨ψ(θ)|H|ψ(θ)⟩ + β|⟨ψ(θ)|ψ₀⟩|²
                                ↑ overlap penalty

The penalty term β|⟨ψ|ψ₀⟩|² adds a large energy cost whenever the new state overlaps with the ground state. This forces the optimizer to find a state that is both low-energy and orthogonal to |ψ₀⟩ — which is exactly the first excited state.

CODECopy
                    Without penalty        With penalty (VQD)
Energy              ┌──────────────┐       ┌──────────────┐
  ↑                 │              │       │              │
  │     E₁ ─────   │  ψ₁ ─────   │       │  ψ₁ ← found │
  │                 │              │       │              │
  │                 │              │       │  ╱╱╱╱╱╱╱╱╱╱  │ ← penalty wall
  │     E₀ ─────   │  ψ₀ ← both  │       │  ψ₀ (known)  │
  │                 │    collapse  │       │              │
  └─────────────    └──────────────┘       └──────────────┘

Step 3: Choosing β

The penalty coefficient β must be larger than the energy gap (E₁ - E₀). If β is too small, the optimizer finds it cheaper to accept the overlap penalty than to move to the excited state — it collapses back to the ground state.

The Ansatz

We use a 2-layer UCCSD-inspired ansatz with RY, RZ, and CX gates:

CODECopy
Layer 1:
     ┌──────────┐┌──────────┐
q_0: ┤ RY(θ₀)  ├┤ RZ(θ₁)  ├──■──
     ├──────────┤├──────────┤┌─┴─┐
q_1: ┤ RY(θ₂)  ├┤ RZ(θ₃)  ├┤ X ├
     └──────────┘└──────────┘└───┘

Layer 2:
     ┌──────────┐┌──────────┐
q_0: ┤ RY(θ₄)  ├┤ RZ(θ₅)  ├──■──
     ├──────────┤├──────────┤┌─┴─┐
q_1: ┤ RY(θ₆)  ├┤ RZ(θ₇)  ├┤ X ├
     └──────────┘└──────────┘└───┘

The RZ gates are critical for excited states. Pure RY ansatze can represent only real-amplitude states, but excited state wavefunctions can have complex amplitudes in the computational basis. The RZ gates add the phase control needed to reach these states.

The Math

VQD Cost Function

For the k-th excited state:

L_k(θ) = ⟨ψ(θ)|H|ψ(θ)⟩ + Σᵢ₌₀ᵏ⁻¹ βᵢ |⟨ψ(θ)|ψᵢ⟩|²

The overlap |⟨ψ(θ)|ψᵢ⟩|² is computed using the statevector inner product in simulation:

overlap = |⟨ψ₁|ψ₂⟩|² = |Σⱼ α*ⱼ βⱼ|²

On real hardware, you would use a swap test — an ancilla-controlled circuit that maps the overlap to a measurement probability:

CODECopy
                ┌───┐           ┌───┐
ancilla: |0⟩ ──┤ H ├──■────────┤ H ├── Measure
                └───┘  │        └───┘
                     ┌─┴──┐
  |ψ₁⟩: ────────────┤SWAP├─────────────
  |ψ₂⟩: ────────────┤    ├─────────────
                     └────┘

P(ancilla = 0) = (1 + |⟨ψ₁|ψ₂⟩|²) / 2

H2 Hamiltonian

The full 2-qubit Hamiltonian (including Y₀Y₁ term) is:

H = -1.0523·I + 0.3979·Z₀ - 0.3979·Z₁ - 0.0112·Z₀Z₁ + 0.1809·X₀X₁ + 0.1809·Y₀Y₁

Exact eigenvalues (electronic only):

• E₀ = -1.9153 Ha (bonding: superposition of |01⟩ and |10⟩)
• E₁ = E₂ = -1.0635 Ha (degenerate pair: |00⟩ and |11⟩ components)
• E₃ = -0.1669 Ha (antibonding: orthogonal superposition of |01⟩ and |10⟩)

Adding nuclear repulsion (V_nn = 0.7199 Ha):

• Ground state total: -1.1954 Ha
• First excited total: -0.3436 Ha
• Gap: 0.8518 Ha ≈ 23.2 eV (ultraviolet)

Expected Output

Running the Circuit

PYTHONCopy
from circuit import run_circuit, optimize_vqd, verify_vqd

# Quick single evaluation
result = run_circuit()
print(f"Energy: {result['energy']:.4f} Ha")

# Full VQD: ground + excited state
opt = optimize_vqd(max_iterations=100, shots=2048, seed=42)
print(f"Ground:  {opt['ground_state']['energy']:.4f} Ha")
print(f"Excited: {opt['excited_state']['energy']:.4f} Ha")
print(f"Gap:     {opt['gap']['computed']:.4f} Ha")
print(f"Overlap: {opt['overlap']:.4f}")

# Verification suite
v = verify_vqd()
for check in v["checks"]:
    status = "PASS" if check["passed"] else "FAIL"
    print(f"[{status}] {check['name']}: {check['detail']}")

Try It Yourself

1. Break β: Set beta=0.5 (below the 0.85 Ha gap) and run optimize_vqd(). The excited state should collapse back to the ground state — the overlap penalty is too weak.

2. Find the second excited state: After finding |ψ₀⟩ and |ψ₁⟩, run a third optimization with penalties against both. Use excited_state_cost() with lower_state_thetas=[ground_theta, excited_theta] and betas=[3.0, 3.0].

3. RY-only ansatz: Remove the RZ gates from the ansatz (set RZ angles to 0). Can it still find the excited state? Why or why not?

4. Swap test overlap: Replace the statevector overlap in compute_overlap() with a swap test circuit. How many shots do you need for the overlap estimate to be accurate within ±0.05?

5. Bond stretching: Change the Hamiltonian coefficients to R = 2.0 A (stretched bond). The gap shrinks — does β = 3.0 still work, or is it now too large?

What's Next

• LiH Excited States — Scale VQD to 4 qubits with lithium hydride
• Multistate VQE — Find multiple states simultaneously (no sequential penalty)
• SSVQE — Subspace-Search VQE: alternative to VQD using orthogonal initial states
• VQD Algorithm — Generalized VQD for arbitrary k-th excited state

Comparison with Other Excited State Methods

Applications

References

• Higgott, O., Wang, D., Brierley, S. (2019). "Variational Quantum Computation of Excited States." Quantum 3, 156. DOI: 10.22331/q-2019-07-01-156
• Peruzzo, A. et al. (2014). "A variational eigenvalue solver on a photonic quantum processor." Nature Communications 5, 4213. DOI: 10.1038/ncomms5213
• McClean, J.R. et al. (2016). "The theory of variational hybrid quantum-classical algorithms." New Journal of Physics 18, 023023. DOI: 10.1088/1367-2630/18/2/023023
• Nakanishi, K.M., Mitarai, K., Fujii, K. (2019). "Subspace-search variational quantum eigensolver for excited states." Physical Review Research 1, 033062. DOI: 10.1103/PhysRevResearch.1.033062