ARPA-E Quantum Simulation Program for Fusion Plasma Physics

The Problem

Fusion energy (the same process that powers the sun) offers effectively unlimited clean energy if a burning plasma can be sustained in a magnetic confinement device (tokamak or stellarator) long enough to produce more energy than is consumed to heat it. The primary obstacle to achieving net energy gain (Q > 1) in current devices is turbulent transport: microscale plasma instabilities driven by temperature and density gradients generate turbulent eddies that rapidly transport energy from the hot plasma core to the cooler edge, cooling the plasma faster than fusion heating can sustain it.

Classical simulation of plasma turbulence uses the gyrokinetic model: a 5-dimensional phase-space distribution function (3 spatial + 2 velocity dimensions) for each ion and electron species, evolved under the gyrokinetic Vlasov-Maxwell equations. A single turbulence simulation for a realistic tokamak geometry requires petascale computing (50,000+ CPU cores, weeks of runtime on Frontier or Summit) to resolve the relevant spatial and temporal scales. This computational cost limits the number of plasma configurations and heating scenarios that can be explored for ITER and future power plant designs. Quantum simulation of the gyrokinetic equation could, in principle, achieve an exponential speedup by encoding the plasma distribution function in a quantum state, a 2^n dimensional Hilbert space accessible with n qubits.

Gyrokinetic Hamiltonian Encoding

The gyrokinetic equation in the delta-f formulation tracks the perturbation delta-f to a Maxwellian background distribution. In a shear-slab geometry with Fourier decomposition in the perpendicular directions, the plasma dynamics reduce to a set of coupled oscillators labeled by wave vector k. The Hamiltonian for each k-mode pair involves the ion acoustic wave, drift wave, and ion temperature gradient (ITG) instability terms.

For a minimal system demonstrating the essential physics (2 species: ions and electrons, and 4 Fourier modes) the Hilbert space has dimension 2^8 = 256, requiring 8 qubits. The Hamiltonian is expressed as a sum of Pauli operators via the Jordan-Wigner transformation of the second-quantized gyrokinetic field operators.

import numpy as np
import pennylane as qml
from pennylane import numpy as pnp

# Gyrokinetic Hamiltonian for slab geometry
# 2 species (ions, electrons), 4 Fourier modes each -> 8 qubits
# H = H_drift + H_ITG + H_coupling

n_qubits = 8  # 4 modes x 2 species
dev = qml.device("default.qubit", wires=n_qubits)

def gyrokinetic_hamiltonian(
    v_drift_i, v_drift_e, eta_i, k_perp_values, coupling_strength
):
    """
    Build gyrokinetic Hamiltonian as PennyLane Hamiltonian.
    v_drift_i/e: ion/electron diamagnetic drift velocities (m/s)
    eta_i: ion temperature gradient parameter (eta = L_n / L_T)
    k_perp_values: perpendicular wavenumbers for each mode (m^-1)
    coupling_strength: ion-electron coupling coefficient
    """
    coeffs = []
    obs = []

    # Single-mode terms: drift wave dispersion for each mode
    for mode in range(4):
        k = k_perp_values[mode]
        # Ion drift wave frequency ~ k * v_drift_i * (1 + eta_i * k^2 * rho_i^2)
        rho_i = 0.01  # ion Larmor radius in meters
        omega_mode_i = k * v_drift_i * (1 + eta_i * k**2 * rho_i**2)
        omega_mode_e = k * v_drift_e

        # Encode as Z operators on ion qubit (modes 0-3) and electron qubit (modes 4-7)
        coeffs.append(omega_mode_i)
        obs.append(qml.PauliZ(mode))         # ion mode

        coeffs.append(omega_mode_e)
        obs.append(qml.PauliZ(mode + 4))     # electron mode

    # Coupling terms: ion-electron interaction via XX coupling
    for mode in range(4):
        k = k_perp_values[mode]
        g = coupling_strength * k
        coeffs.append(g)
        obs.append(qml.PauliX(mode) @ qml.PauliX(mode + 4))

        coeffs.append(g)
        obs.append(qml.PauliY(mode) @ qml.PauliY(mode + 4))

    return qml.Hamiltonian(coeffs, obs)

# ITER-like parameters (normalized)
k_perp = np.array([0.1, 0.2, 0.3, 0.4])  # k_perp * rho_i values
H_gyrokinetic = gyrokinetic_hamiltonian(
    v_drift_i=1e4,   # ion diamagnetic drift (m/s)
    v_drift_e=3e4,   # electron drift
    eta_i=2.0,       # ITG instability threshold: eta_i > 2/3
    k_perp_values=k_perp,
    coupling_strength=500.0
)
print(f"Gyrokinetic Hamiltonian terms: {len(H_gyrokinetic.ops)}")
print(f"Number of qubits required: {n_qubits}")

ADAPT-VQE for Plasma Ground State

The plasma ground state (lowest energy configuration of the turbulent fluctuations) determines the nonlinear saturation amplitude of the turbulence, which in turn sets the turbulent transport level. ADAPT-VQE (Adaptive Derivative-Assembled Pseudo-Trotter Variational Quantum Eigensolver) builds the ansatz circuit adaptively by selecting the operator at each step that maximally reduces the energy gradient, producing shallower circuits than fixed UCCSD ansatze for systems with sparse Hamiltonians.

ARPA-E-funded teams at MIT, Princeton, and the University of Colorado applied ADAPT-VQE to the 8-qubit gyrokinetic Hamiltonian on both IBM Eagle (via cloud access) and IonQ Aria (trapped ion, lower gate error rates for deep circuits).

import pennylane as qml
from pennylane import numpy as pnp
import numpy as np

dev = qml.device("default.qubit", wires=n_qubits)

# ADAPT-VQE operator pool: single and double excitations between modes
def build_operator_pool(n_qubits):
    """Generate operator pool for ADAPT-VQE on gyrokinetic system."""
    pool = []
    # Single-qubit rotations
    for i in range(n_qubits):
        for pauli in [qml.PauliX, qml.PauliY]:
            pool.append([(1.0, pauli(i))])

    # Two-qubit entangling operators: XX and YY between ion-electron pairs
    for mode in range(4):
        pool.append([
            (0.5, qml.PauliX(mode) @ qml.PauliX(mode + 4)),
            (-0.5, qml.PauliY(mode) @ qml.PauliY(mode + 4))
        ])
        pool.append([
            (0.5, qml.PauliY(mode) @ qml.PauliX(mode + 4)),
            (0.5, qml.PauliX(mode) @ qml.PauliY(mode + 4))
        ])
    return pool

def adapt_vqe_step(H, current_state_params, operator_pool, threshold=1e-3):
    """
    One ADAPT-VQE iteration: select operator with largest gradient, add to ansatz.
    Returns selected operator index and gradient magnitude.
    """
    gradients = []
    for op_terms in operator_pool:
        # Compute gradient: <psi| [H, A_mu] |psi> for each pool operator A_mu
        grad_magnitude = np.random.exponential(0.1)  # placeholder for actual QPU measurement
        gradients.append(grad_magnitude)

    best_idx = np.argmax(gradients)
    best_gradient = gradients[best_idx]
    converged = best_gradient < threshold

    return best_idx, best_gradient, converged

operator_pool = build_operator_pool(n_qubits)
print(f"Operator pool size: {len(operator_pool)}")

# ADAPT-VQE main loop
params = []
selected_operators = []
n_adapt_iterations = 10

for iteration in range(n_adapt_iterations):
    op_idx, grad, converged = adapt_vqe_step(
        H_gyrokinetic, params, operator_pool
    )
    selected_operators.append(op_idx)
    params.append(0.01)  # initialize new parameter near zero
    print(f"Iteration {iteration+1}: selected operator {op_idx}, gradient={grad:.4f}")
    if converged:
        print("ADAPT-VQE converged")
        break

print(f"Ansatz depth (operators selected): {len(selected_operators)}")

Quantum Dynamics for Turbulent Transport

Beyond ground state energy, the transport coefficient (thermal diffusivity chi) is determined by the time evolution of the turbulent fluctuations. ARPA-E teams used quantum simulation of the time-dependent Schrodinger equation for the gyrokinetic Hamiltonian, implementing Suzuki-Trotter decomposition to second order. The transport coefficient is extracted from the long-time correlation function of the turbulent heat flux operator Q = sum_k k * delta-f_k * phi_k, where phi_k is the electrostatic potential fluctuation.

For the 8-qubit, 2-species, 4-mode system, the Trotter simulation ran for 500 time steps with a step size of dt = 0.01 (normalized to the ion cyclotron period). The computed thermal diffusivity was chi_i = 0.4 chi_GB (in gyro-Bohm units), consistent with the linear gyrokinetic prediction for this parameter regime. Full nonlinear gyrokinetic simulation of an ITER plasma cross-section requires encoding approximately 10^6 Fourier modes per species across a 2D flux surface geometry, corresponding to millions of logical qubits. Current hardware supports only the smallest proof-of-concept problems, but the algorithmic framework for quantum gyrokinetic simulation is now established and demonstrated. ARPA-E’s 8 funded teams are publishing the theoretical foundations and small-scale demonstrations that will guide hardware roadmaps over the next decade.

Learn more: Qiskit Reference