Analog Hamiltonian Simulation with QuEra Aquila

What Is Analog Hamiltonian Simulation?

Gate-based quantum computing decomposes an algorithm into a sequence of discrete unitary operations (gates) applied to individual qubits. Analog Hamiltonian simulation takes a fundamentally different approach: it evolves a physical system under a continuously driven Hamiltonian that you shape by tuning external control fields in real time. There is no gate decomposition, no circuit compilation, and no concept of “circuit depth.” Instead, you specify (1) where the atoms sit and (2) how the laser fields change over time, and the physics of the system does the rest.

QuEra’s Aquila processor, available through Amazon Braket, uses arrays of neutral rubidium-87 atoms held in optical tweezers. Each atom has two relevant states: the electronic ground state |g> and a highly excited Rydberg state |r>. A global laser field couples these two states. When two atoms are close enough together, the strong van der Waals interaction between their Rydberg states prevents both atoms from being excited simultaneously. This is the Rydberg blockade, and it is the core mechanism that generates entanglement and many-body correlations in the system.

Analog Hamiltonian simulation on Aquila is well suited for three classes of problems:

1. Quantum phase transitions: sweeping control parameters across a critical point to observe ordered phases of matter.
2. Combinatorial optimization: encoding constraint satisfaction problems (like Maximum Independent Set) into the blockade geometry.
3. Many-body physics: studying dynamics that are intractable on classical computers, such as thermalization and quantum scarring.

This tutorial walks through the physics, the API, and the practical details of programming Aquila through Braket. By the end, you will be able to design atom arrangements, build adiabatic pulse sequences, interpret measurement results, compute spatial correlations, benchmark against classical solvers, and submit jobs to real hardware.

The Rydberg Hamiltonian

Aquila implements the following Hamiltonian:

H = (Omega/2) * sum_i (e^(i*phi)|g_i><r_i| + h.c.) - Delta * sum_i n_i + sum_{i<j} V_ij * n_i * n_j

Each term has a distinct physical role:

• (Omega/2) * sum_i (e^(i*phi)|g_i><r_i| + h.c.): The driving term. The global laser with Rabi frequency Omega and phase phi coherently couples the ground state |g> to the Rydberg state |r> on every atom. This term creates superpositions of |g> and |r>.

• -Delta * sum_i n_i: The detuning term. Here n_i = |r_i><r_i| is the number operator that counts whether atom i is in the Rydberg state. When Delta is large and negative, atoms prefer to stay in |g> (the detuning penalty for being in |r> is positive). When Delta is large and positive, atoms prefer to be in |r> (they gain energy by being excited). The detuning acts as a chemical potential for Rydberg excitations.

• sum_{i<j} V_ij * n_i * n_j: The interaction term. V_ij = C6 / |r_i - r_j|^6 is the van der Waals interaction between atoms i and j. This interaction is repulsive (V_ij > 0 for the Rydberg states used in Aquila), so it penalizes configurations where two nearby atoms are both in |r>. This is the origin of the Rydberg blockade.

The key insight: you control Omega(t), Delta(t), and phi(t) as functions of time, but the interaction V_ij is fixed by the atom positions. By choosing atom positions wisely, you encode the problem structure into the interaction graph. By choosing the pulse shapes wisely, you steer the system toward the ground state of the final Hamiltonian.

Rydberg Blockade Physics in Depth

The van der Waals interaction V_ij = C6 / r^6 drops off steeply with distance. For rubidium-87 atoms in the 70S Rydberg state, the C6 coefficient is approximately 862,690 atomic units. Converting to SI: 1 atomic unit of C6 = E_h * a_0^6 where E_h = 4.36 x 10^-18 J and a_0 = 5.29 x 10^-11 m. In practice, C6 ≈ 5.42 x 10^-24 Jm^6 (equivalently, about 2pi * 862,690 MHz * micrometer^6 in the frequency-distance units commonly used in the field).

The blockade radius r_b is defined as the distance at which the interaction energy equals the driving strength:

V(r_b) = C6 / r_b^6 = Omega / 2

Solving for r_b:

r_b = (2 * C6 / Omega)^(1/6)

Let’s compute this for a typical Rabi frequency:

import numpy as np

# C6 coefficient for Rb-87, 70S Rydberg state, in frequency*distance^6 units
# C6 ≈ 2*pi * 862690 MHz * micrometer^6
C6 = 2 * np.pi * 862690e6 * (1e-6)**6  # Convert to rad/s * m^6

# Rabi frequency: Omega = 2*pi * 2.5 MHz
omega = 2 * np.pi * 2.5e6  # rad/s

# Blockade radius
r_b = (2 * C6 / omega) ** (1/6)
print(f"Blockade radius: {r_b * 1e6:.2f} micrometers")

# Interaction at 5.5 micrometer spacing (our lattice spacing)
r_nn = 5.5e-6  # meters
V_nn = C6 / r_nn**6
print(f"Nearest-neighbor interaction V_nn = 2*pi * {V_nn / (2 * np.pi * 1e6):.2f} MHz")
print(f"V_nn / Omega = {V_nn / omega:.2f}")
print(f"Nearest neighbors ARE within blockade radius: {r_nn < r_b}")

# Interaction at diagonal distance (5.5*sqrt(2) micrometers)
r_diag = 5.5e-6 * np.sqrt(2)
V_diag = C6 / r_diag**6
print(f"\nDiagonal interaction V_diag = 2*pi * {V_diag / (2 * np.pi * 1e6):.2f} MHz")
print(f"V_diag / Omega = {V_diag / omega:.2f}")
print(f"Diagonal neighbors within blockade radius: {r_diag < r_b}")

For Omega = 2*pi * 2.5 MHz, the blockade radius is approximately 8.5 micrometers. At our chosen lattice spacing of 5.5 micrometers, nearest-neighbor atoms have V_nn much larger than Omega, so they are strongly blockaded. Even diagonal neighbors at 5.5 * sqrt(2) ≈ 7.8 micrometers are within the blockade radius. This means that on a 3x3 grid with 5.5 micrometer spacing, no two adjacent atoms (horizontally, vertically, or diagonally) can both be in |r>. The blockade graph is the king graph (8-connected grid), not just the 4-connected grid.

If you want only nearest neighbors to be blockaded (not diagonals), increase the lattice spacing to around 7 micrometers. Alternatively, reduce Omega to shrink r_b. The interplay between lattice geometry and blockade radius determines the effective interaction graph, which in turn determines what optimization problem the system solves.

The Adiabatic Sweep: Intuition and Theory

The adiabatic theorem from quantum mechanics states: if you change the Hamiltonian slowly enough compared to the energy gap between the ground state and first excited state, the system remains in the instantaneous ground state throughout the evolution.

For the Rydberg Hamiltonian, the adiabatic sweep works as follows:

1. Initial state (Delta << -Omega): When the detuning is large and negative, the energy cost of exciting any atom to |r> far exceeds the driving strength. The ground state is |ggg…g>, all atoms in the ground state. This state is easy to prepare: it is the natural starting point when no laser is applied.

2. Intermediate regime (|Delta| ~ Omega): The driving competes with the detuning. Atoms begin to transition between |g> and |r>. The system passes through a quantum phase transition, where quantum fluctuations are maximal and the energy gap between the ground and first excited state is smallest. This is the bottleneck for adiabaticity.

3. Final state (Delta >> +Omega): The detuning strongly favors Rydberg excitation. Every atom “wants” to be in |r>, but the blockade interaction prevents neighbors from both being excited. The ground state is the configuration that maximizes the number of Rydberg excitations subject to the blockade constraint. This is exactly the Maximum Independent Set of the blockade graph.

The quality of the adiabatic preparation depends on how slowly you sweep through the critical point. The adiabatic condition requires:

|d(Delta)/dt| << gap^2 / |<excited| dH/dDelta |ground>|

where gap is the minimum energy gap during the evolution. For many-body systems, this gap typically shrinks polynomially (or even exponentially for first-order transitions) with system size. For the small 3x3 grids in this tutorial, a 4-microsecond sweep is generally sufficient.

The pulse protocol consists of three stages:

1. Ramp up Omega while Delta stays negative. This gradually introduces quantum fluctuations without exciting atoms (the large negative detuning keeps them in |g>).
2. Sweep Delta from negative to positive while Omega is at its maximum. This drives the system through the phase transition.
3. Ramp down Omega while Delta stays positive. This projects the quantum state onto the classical basis of Rydberg/ground configurations. Without this ramp-down, the measurement basis would not align with the independent set encoding.

Maximum Independent Set: The Optimization Connection

The connection between Rydberg blockade physics and the Maximum Independent Set (MIS) problem is the primary reason Aquila is interesting for optimization.

Definition: Given a graph G = (V, E), an independent set is a subset S of vertices such that no two vertices in S share an edge. The Maximum Independent Set is the largest such subset.

MIS is NP-hard in general. No known classical algorithm solves it in polynomial time for arbitrary graphs.

The Rydberg encoding: Construct the blockade graph by placing atoms at specific positions and choosing Omega such that:

• There is an edge between atoms i and j if and only if their distance r_ij < r_b (the blockade radius).
• The ground state of the Rydberg Hamiltonian at large positive Delta is the MIS of this blockade graph.

For the 3x3 grid with 5.5 micrometer spacing and Omega = 2*pi * 2.5 MHz, the blockade graph connects all atoms within about 8.5 micrometers of each other. This includes nearest neighbors (5.5 micrometers) and diagonal neighbors (7.8 micrometers), producing the king graph on a 3x3 grid.

The MIS of the 3x3 king graph has size 3 (you can verify: place excitations on alternating atoms such that no two share a horizontal, vertical, or diagonal edge). However, if you increase the spacing to about 7 micrometers so that only nearest neighbors (not diagonals) are blockaded, the blockade graph becomes the 4-connected grid, and the MIS has size 5 (the checkerboard pattern).

# Visualize the two possible blockade graphs for a 3x3 grid

import itertools

def build_blockade_graph(positions, blockade_radius):
    """Build adjacency list for the blockade graph."""
    n = len(positions)
    adj = {i: [] for i in range(n)}
    for i in range(n):
        for j in range(i + 1, n):
            dx = positions[i][0] - positions[j][0]
            dy = positions[i][1] - positions[j][1]
            dist = np.sqrt(dx**2 + dy**2)
            if dist < blockade_radius:
                adj[i].append(j)
                adj[j].append(i)
    return adj

def find_mis_brute_force(adj, n):
    """Find maximum independent set by brute force (only for small graphs)."""
    best_set = []
    for size in range(n, 0, -1):
        for subset in itertools.combinations(range(n), size):
            valid = True
            for v in subset:
                for u in subset:
                    if u != v and u in adj[v]:
                        valid = False
                        break
                if not valid:
                    break
            if valid:
                return list(subset)
    return best_set

# 3x3 grid positions in micrometers
spacing_tight = 5.5e-6
positions_3x3 = []
for row in range(3):
    for col in range(3):
        positions_3x3.append((col * spacing_tight, row * spacing_tight))

# Blockade radius for Omega = 2*pi*2.5 MHz
r_b = 8.5e-6

# King graph (diagonals included)
adj_king = build_blockade_graph(positions_3x3, r_b)
mis_king = find_mis_brute_force(adj_king, 9)
print("King graph (5.5 um spacing, r_b=8.5 um):")
print(f"  Edges per vertex: {[len(adj_king[i]) for i in range(9)]}")
print(f"  MIS size: {len(mis_king)}, vertices: {mis_king}")

# 4-connected graph (increase spacing so diagonals are outside blockade)
spacing_wide = 7.0e-6
positions_wide = []
for row in range(3):
    for col in range(3):
        positions_wide.append((col * spacing_wide, row * spacing_wide))

# Reduce blockade radius (or equivalently, increase spacing)
# Diagonal distance = 7*sqrt(2) = 9.9 um, which is > r_b = 8.5 um
adj_grid = build_blockade_graph(positions_wide, r_b)
mis_grid = find_mis_brute_force(adj_grid, 9)
print(f"\n4-connected grid (7.0 um spacing, r_b=8.5 um):")
print(f"  Edges per vertex: {[len(adj_grid[i]) for i in range(9)]}")
print(f"  MIS size: {len(mis_grid)}, vertices: {mis_grid}")

# Print the MIS as a grid
def print_grid(mis_vertices, cols=3, rows=3):
    for r in range(rows):
        row_str = ""
        for c in range(cols):
            idx = r * cols + c
            row_str += "R " if idx in mis_vertices else "g "
        print(f"  {row_str}")

print("\nMIS pattern (king graph):")
print_grid(mis_king)
print("\nMIS pattern (4-connected grid, checkerboard):")
print_grid(mis_grid)

The choice of lattice spacing relative to the blockade radius directly determines the optimization problem that the hardware solves. This is a powerful design knob.

Setting Up Atom Positions

from braket.ahs import AtomArrangement, DrivingField, AnalogHamiltonianSimulation
from braket.ahs import TimeSeries
from braket.devices import LocalSimulator
import numpy as np

# Create a 3x3 square lattice with spacing 5.5 micrometers
# This spacing is chosen so nearest neighbors are within the blockade radius
lattice_spacing = 5.5e-6  # meters

register = AtomArrangement()
for row in range(3):
    for col in range(3):
        x = col * lattice_spacing
        y = row * lattice_spacing
        register.add([x, y])

print(f"Number of atoms: {len(register)}")

Atom positions are specified in SI units (meters). Aquila supports up to 256 atom sites in custom 2D arrangements, though not all sites need to be filled.

Aquila Hardware Specifications and Constraints

Before designing your pulse sequence, you need to know what the hardware allows. The Aquila device has the following specifications and constraints:

Parameter	Value
Maximum atom sites	256
Minimum inter-atom spacing	4 micrometers
Field of view	75 x 76 micrometers
Atom loading probability	> 99% per site
Single-atom state detection fidelity	> 99.5%
Coherence time (T2)	1 to 2 microseconds
Maximum Rabi frequency (Omega)	2*pi * 4 MHz
Maximum detuning magnitude (|Delta|)	2*pi * 20 MHz
Minimum total evolution time	0.4 microseconds
Maximum total evolution time	4 microseconds
Atomic species	Rb-87
Rydberg state	70S

These specifications constrain your experiment design in several ways:

• Atom arrangement: With a 75 x 76 micrometer field of view and 4 micrometer minimum spacing, you could in principle fit a grid of about 19 x 19 = 361 sites. The 256-site limit is more restrictive, but still allows large and complex geometries.
• Pulse timing: The 4-microsecond maximum evolution time limits how adiabatic your sweep can be. For large systems where the energy gap shrinks, this ceiling can be the bottleneck for solution quality.
• Drive amplitude: The maximum Omega of 2*pi * 4 MHz limits the blockade radius. For larger blockade radii, you must use lower Omega, which further constrains the adiabatic sweep speed.
• Coherence: The T2 coherence time of 1 to 2 microseconds means that for evolution times approaching or exceeding T2, decoherence will degrade the fidelity of your prepared state. This is a real constraint for 4-microsecond sweeps.

Pulse Shape Design and API Constraints

The piecewise-linear TimeSeries is the standard control scheme in the Braket AHS API. The API enforces several constraints on the pulse shapes. Here is a reference implementation with each constraint annotated:

import numpy as np
from braket.ahs import TimeSeries, DrivingField

# Hardware limits (from Aquila device specifications)
OMEGA_MAX = 2 * np.pi * 4e6       # Maximum Rabi frequency: 2*pi * 4 MHz
DELTA_MAX = 2 * np.pi * 20e6      # Maximum |Delta|: 2*pi * 20 MHz
T_MIN = 0.4e-6                     # Minimum total time: 0.4 microseconds
T_MAX = 4.0e-6                     # Maximum total time: 4.0 microseconds
MIN_SPACING = 4.0e-6               # Minimum inter-atom distance: 4 micrometers

# Chosen parameters (within hardware limits)
t_total = 4e-6                     # Use the full 4 microseconds for best adiabaticity
omega_drive = 2 * np.pi * 2.5e6   # 2.5 MHz, well under the 4 MHz max
delta_start = -2 * np.pi * 5e6    # -5 MHz, within the +/-20 MHz range
delta_end = 2 * np.pi * 5e6       # +5 MHz

# Constraint: Omega MUST start at 0 and end at 0
# The hardware cannot instantaneously turn on or off the laser.
# A discontinuity at t=0 or t=T would require infinite bandwidth.
t_ramp = 0.5e-6  # 500 ns ramp time

amplitude = TimeSeries()
amplitude.put(0.0, 0.0)                    # Start at zero (required)
amplitude.put(t_ramp, omega_drive)          # Ramp up over 500 ns
amplitude.put(t_total - t_ramp, omega_drive)  # Hold at constant Omega
amplitude.put(t_total, 0.0)                 # Ramp back to zero (required)

# Detuning: linear sweep from negative to positive
# No requirement that Delta starts or ends at zero
detuning = TimeSeries()
detuning.put(0.0, delta_start)
detuning.put(t_total, delta_end)

# Phase: constant at zero for standard protocols
# Phase control can be used for more advanced pulse sequences
phase = TimeSeries()
phase.put(0.0, 0.0)
phase.put(t_total, 0.0)

drive = DrivingField(amplitude=amplitude, detuning=detuning, phase=phase)

Real experiments at QuEra often use smoother pulse shapes (Gaussian ramps, sinusoidal profiles) to minimize leakage to non-target states. The piecewise-linear TimeSeries can approximate smooth shapes by using many short segments. For example, a Gaussian ramp-up:

def gaussian_ramp(t_start, t_end, n_points, amplitude_max):
    """Create a Gaussian ramp-up using piecewise-linear segments."""
    ts = TimeSeries()
    times = np.linspace(t_start, t_end, n_points)
    sigma = (t_end - t_start) / 4  # 4-sigma width
    center = t_start
    for t in times:
        # Gaussian CDF-like profile from 0 to amplitude_max
        val = amplitude_max * (1 - np.exp(-0.5 * ((t - t_start) / sigma) ** 2))
        ts.put(float(t), float(val))
    return ts

Building and Running the Program

from braket.ahs import AtomArrangement, DrivingField, AnalogHamiltonianSimulation
from braket.ahs import TimeSeries
from braket.devices import LocalSimulator
import numpy as np

# Set up the 3x3 lattice
lattice_spacing = 5.5e-6
register = AtomArrangement()
for row in range(3):
    for col in range(3):
        register.add([col * lattice_spacing, row * lattice_spacing])

# Build the drive (same parameters as above)
t_total = 4e-6
omega_drive = 2 * np.pi * 2.5e6
delta_start = -2 * np.pi * 5e6
delta_end = 2 * np.pi * 5e6
t_ramp = 0.5e-6

amplitude = TimeSeries()
amplitude.put(0.0, 0.0)
amplitude.put(t_ramp, omega_drive)
amplitude.put(t_total - t_ramp, omega_drive)
amplitude.put(t_total, 0.0)

detuning = TimeSeries()
detuning.put(0.0, delta_start)
detuning.put(t_total, delta_end)

phase = TimeSeries()
phase.put(0.0, 0.0)
phase.put(t_total, 0.0)

drive = DrivingField(amplitude=amplitude, detuning=detuning, phase=phase)

# Construct the AHS program
ahs_program = AnalogHamiltonianSimulation(register=register, hamiltonian=drive)

# Run on the local simulator
local_sim = LocalSimulator("braket_ahs")
result = local_sim.run(ahs_program, shots=1000).result()

# Each shot produces a measurement with pre_sequence and post_sequence
measurements = result.measurements

# Count Rydberg excitations per shot
# IMPORTANT: post_sequence uses 0 for Rydberg-excited and 1 for ground state
# This is counterintuitive. 0 means "not detected in the ground state",
# which indicates the atom is in the Rydberg state (or was lost).
rydberg_counts = []
for shot in measurements:
    n_rydberg = sum(1 for s in shot.post_sequence if s == 0)
    rydberg_counts.append(n_rydberg)

avg_rydberg = np.mean(rydberg_counts)
print(f"Average Rydberg excitations per shot: {avg_rydberg:.2f} out of 9 atoms")
print(f"Standard deviation: {np.std(rydberg_counts):.2f}")

# Histogram of Rydberg counts
from collections import Counter
count_hist = Counter(rydberg_counts)
for n in sorted(count_hist.keys()):
    frac = count_hist[n] / len(rydberg_counts)
    bar = "#" * int(frac * 50)
    print(f"  {n} Rydberg atoms: {count_hist[n]:4d} shots ({frac:.1%}) {bar}")

Interpreting Results: Phase Transition Signature

When the detuning sweeps from far below resonance to far above, the system transitions from the ground state (all atoms in |g>) to a Rydberg-ordered phase. Due to the blockade constraint, neighboring atoms cannot both be in |r>, so the system forms patterns that maximize Rydberg excitations while respecting the blockade.

def check_blockade_violations(measurements, register, blockade_radius=7e-6):
    """Count shots with nearest-neighbor double Rydberg excitations."""
    positions = [(site[0], site[1]) for site in register]
    violations = 0

    for shot in measurements:
        rydberg_sites = [i for i, s in enumerate(shot.post_sequence) if s == 0]
        shot_violated = False
        for idx_a in range(len(rydberg_sites)):
            for idx_b in range(idx_a + 1, len(rydberg_sites)):
                a, b = rydberg_sites[idx_a], rydberg_sites[idx_b]
                dist = np.sqrt((positions[a][0] - positions[b][0])**2 +
                               (positions[a][1] - positions[b][1])**2)
                if dist < blockade_radius:
                    shot_violated = True
                    break
            if shot_violated:
                break
        if shot_violated:
            violations += 1

    return violations

n_violations = check_blockade_violations(result.measurements, register)
print(f"Blockade violations: {n_violations} out of {len(result.measurements)} shots")
print(f"Violation rate: {n_violations / len(result.measurements):.1%}")

If the adiabatic sweep is slow enough, you should see very few blockade violations. The ratio of Rydberg-excited atoms and the suppression of nearest-neighbor co-excitations are the hallmarks of the ordered phase.

Spatial Correlation Functions

Counting Rydberg atoms tells you about the average density, but spatial correlations reveal the structure of the ordered phase. The connected correlation function measures how the excitation of one atom is correlated with another at distance r:

C(r) = <n_i * n_j> - <n_i> * <n_j>

where the average is taken over all pairs (i, j) with separation r and over all measurement shots.

• C(r) < 0: Anti-correlation. Atoms at this distance tend not to both be excited. This is the blockade signature for short distances.
• C(r) > 0: Positive correlation. Atoms at this distance tend to be excited together. For a checkerboard ordering on a grid, next-nearest neighbors (diagonal atoms) show positive correlation because they belong to the same sublattice.
• C(r) ≈ 0: No correlation. The excitation of one atom provides no information about the other.

def compute_spatial_correlations(measurements, register):
    """Compute the connected correlation function C(r) from measurement data."""
    positions = [(site[0], site[1]) for site in register]
    n_atoms = len(positions)
    n_shots = len(measurements)

    # Convert measurements to binary arrays: 1 = Rydberg, 0 = ground
    # Remember: post_sequence 0 = Rydberg, 1 = ground
    bitstrings = []
    for shot in measurements:
        bits = [1 if s == 0 else 0 for s in shot.post_sequence]
        bitstrings.append(bits)
    bitstrings = np.array(bitstrings)  # shape: (n_shots, n_atoms)

    # Compute single-site averages <n_i>
    avg_n = np.mean(bitstrings, axis=0)  # shape: (n_atoms,)

    # Compute pair correlations <n_i * n_j> for each pair
    # Group by distance
    from collections import defaultdict
    correlations_by_dist = defaultdict(list)

    for i in range(n_atoms):
        for j in range(i + 1, n_atoms):
            dx = positions[i][0] - positions[j][0]
            dy = positions[i][1] - positions[j][1]
            dist = np.sqrt(dx**2 + dy**2)

            # Round distance to nearest 0.1 micrometer for binning
            dist_binned = round(dist * 1e6, 1)  # in micrometers

            # Connected correlation
            ni_nj_avg = np.mean(bitstrings[:, i] * bitstrings[:, j])
            connected = ni_nj_avg - avg_n[i] * avg_n[j]
            correlations_by_dist[dist_binned].append(connected)

    # Average over pairs at the same distance
    distances = sorted(correlations_by_dist.keys())
    avg_correlations = []
    for d in distances:
        avg_c = np.mean(correlations_by_dist[d])
        avg_correlations.append(avg_c)
        print(f"  r = {d:5.1f} um: C(r) = {avg_c:+.4f}  ({len(correlations_by_dist[d])} pairs)")

    return distances, avg_correlations

print("Spatial correlations:")
distances, correlations = compute_spatial_correlations(result.measurements, register)

For a well-prepared blockade-ordered state, you should see:

• Strong negative correlation at the nearest-neighbor distance (5.5 micrometers): the blockade suppresses co-excitation.
• Weak positive or near-zero correlation at the diagonal distance (7.8 micrometers): these atoms are in the same sublattice and tend to be excited together.
• Near-zero correlation at longer distances: the correlations decay.

Local Detuning for Weighted Optimization Problems

The global driving field applies the same Omega and Delta to every atom. But Aquila also supports local detuning: a spatially varying detuning field that assigns different Delta values to different atoms. This is the mechanism for encoding weighted optimization problems.

The local detuning adds a term to the Hamiltonian:

H_local = -sum_i delta_i(t) * n_i

where delta_i(t) is the local detuning on atom i. In the Braket API, this is specified using a ShiftingField object that combines a spatial Pattern (the relative weight per atom) with a TimeSeries (the global time dependence).

from braket.ahs import ShiftingField, Pattern

# Suppose we have a 3x3 grid and want to encode vertex weights.
# Atoms with higher weight should be more likely to be in the MIS.
# We achieve this by giving them larger local detuning (more incentive to be in |r>).

# Example weights for 9 vertices
vertex_weights = [1.0, 2.0, 1.0,
                  2.0, 3.0, 2.0,
                  1.0, 2.0, 1.0]

# Normalize weights to [0, 1] range (Braket requires pattern values in [0, 1])
max_weight = max(vertex_weights)
normalized_weights = [w / max_weight for w in vertex_weights]

# Create the spatial pattern
local_pattern = Pattern(normalized_weights)

# Time dependence: ramp up the local detuning during the sweep
# The local detuning magnitude is the global_max * pattern_value for each atom
local_detuning_max = 2 * np.pi * 2e6  # 2 MHz maximum local detuning

local_ts = TimeSeries()
local_ts.put(0.0, 0.0)              # Start with no local detuning
local_ts.put(t_ramp, 0.0)           # Keep zero during Omega ramp-up
local_ts.put(t_total - t_ramp, local_detuning_max)  # Ramp up with the sweep
local_ts.put(t_total, local_detuning_max)

local_shift = ShiftingField(magnitude=local_ts, pattern=local_pattern)

# Build the AHS program with local detuning
ahs_program_weighted = AnalogHamiltonianSimulation(
    register=register,
    hamiltonian=drive + local_shift  # Add local detuning to the global drive
)

print("AHS program with local detuning created.")
print(f"Atom with highest weight (index 4, center): local detuning = {local_detuning_max/(2*np.pi*1e6):.1f} MHz")
print(f"Atom with lowest weight (index 0, corner): local detuning = {local_detuning_max * normalized_weights[0]/(2*np.pi*1e6):.1f} MHz")

The local detuning biases the optimization: atoms with higher delta_i are more energetically favorable to excite, so they appear more often in the measured independent sets. This transforms the problem from unweighted MIS to Maximum Weight Independent Set (MWIS), which is the more general and practically relevant optimization problem.

Running Parameter Sweeps

The quality of the adiabatic solution depends on the sweep rate. To find the best annealing schedule, run multiple programs with different total evolution times and compare results.

from braket.ahs import (
    AtomArrangement, DrivingField, AnalogHamiltonianSimulation, TimeSeries
)
from braket.devices import LocalSimulator
import numpy as np

# Fixed parameters
lattice_spacing = 5.5e-6
omega_drive = 2 * np.pi * 2.5e6
delta_start = -2 * np.pi * 5e6
delta_end = 2 * np.pi * 5e6

# Build register once
register = AtomArrangement()
for row in range(3):
    for col in range(3):
        register.add([col * lattice_spacing, row * lattice_spacing])

# Sweep over total evolution times
# Longer times are more adiabatic but more susceptible to decoherence
total_times = [1.0e-6, 2.0e-6, 3.0e-6, 4.0e-6]
n_shots = 500

local_sim = LocalSimulator("braket_ahs")
sweep_results = {}

for t_total in total_times:
    t_ramp = min(0.5e-6, t_total * 0.15)  # Ramp is 15% of total, max 500 ns

    amp = TimeSeries()
    amp.put(0.0, 0.0)
    amp.put(t_ramp, omega_drive)
    amp.put(t_total - t_ramp, omega_drive)
    amp.put(t_total, 0.0)

    det = TimeSeries()
    det.put(0.0, delta_start)
    det.put(t_total, delta_end)

    ph = TimeSeries()
    ph.put(0.0, 0.0)
    ph.put(t_total, 0.0)

    drive = DrivingField(amplitude=amp, detuning=det, phase=ph)
    program = AnalogHamiltonianSimulation(register=register, hamiltonian=drive)
    result = local_sim.run(program, shots=n_shots).result()

    # Analyze results
    rydberg_counts = []
    for shot in result.measurements:
        n_r = sum(1 for s in shot.post_sequence if s == 0)
        rydberg_counts.append(n_r)

    avg_r = np.mean(rydberg_counts)
    max_r = max(rydberg_counts)

    # Count how many shots match the optimal MIS configuration
    # For the king graph on 3x3 grid, MIS size is 3
    optimal_count = sum(1 for c in rydberg_counts if c == 3)
    optimal_frac = optimal_count / n_shots

    sweep_results[t_total] = {
        "avg_rydberg": avg_r,
        "max_rydberg": max_r,
        "optimal_fraction": optimal_frac,
    }

    print(f"t_total = {t_total*1e6:.1f} us: avg Rydberg = {avg_r:.2f}, "
          f"max = {max_r}, optimal fraction = {optimal_frac:.1%}")

# Expect: longer sweep times produce higher optimal fractions
# until decoherence effects dominate (not modeled in LocalSimulator)

On the local simulator (which does not model decoherence), longer sweep times should monotonically improve the solution quality. On real hardware, there is an optimal sweep time that balances adiabaticity against decoherence. Finding this sweet spot is a key part of experimental design.

Defects and Position Disorder

Real optical tweezers have finite positioning accuracy. Atoms are not placed at their target positions perfectly; each atom has a small random displacement. This position disorder changes the interaction strengths V_ij and can degrade the quality of the prepared state.

import numpy as np
from braket.ahs import (
    AtomArrangement, DrivingField, AnalogHamiltonianSimulation, TimeSeries
)
from braket.devices import LocalSimulator

def run_with_disorder(lattice_spacing, disorder_std, n_trials=10, n_shots=200):
    """Run AHS simulations with Gaussian position disorder."""
    local_sim = LocalSimulator("braket_ahs")
    t_total = 4e-6
    omega_drive = 2 * np.pi * 2.5e6
    t_ramp = 0.5e-6

    trial_results = []

    for trial in range(n_trials):
        # Create register with noisy positions
        register = AtomArrangement()
        for row in range(3):
            for col in range(3):
                x = col * lattice_spacing + np.random.normal(0, disorder_std)
                y = row * lattice_spacing + np.random.normal(0, disorder_std)
                # Ensure minimum spacing is respected (Aquila requires >= 4 um)
                register.add([x, y])

        # Same drive for all trials
        amp = TimeSeries()
        amp.put(0.0, 0.0)
        amp.put(t_ramp, omega_drive)
        amp.put(t_total - t_ramp, omega_drive)
        amp.put(t_total, 0.0)

        det = TimeSeries()
        det.put(0.0, -2 * np.pi * 5e6)
        det.put(t_total, 2 * np.pi * 5e6)

        ph = TimeSeries()
        ph.put(0.0, 0.0)
        ph.put(t_total, 0.0)

        drive = DrivingField(amplitude=amp, detuning=det, phase=ph)
        program = AnalogHamiltonianSimulation(register=register, hamiltonian=drive)

        result = local_sim.run(program, shots=n_shots).result()

        rydberg_counts = [
            sum(1 for s in shot.post_sequence if s == 0)
            for shot in result.measurements
        ]
        trial_results.append(np.mean(rydberg_counts))

    return trial_results

# Compare different disorder levels
disorder_levels = [0.0, 0.1e-6, 0.3e-6, 0.5e-6]  # in meters

print("Effect of position disorder on Rydberg excitation density:")
for sigma in disorder_levels:
    results = run_with_disorder(5.5e-6, sigma, n_trials=10, n_shots=200)
    mean_r = np.mean(results)
    std_r = np.std(results)
    print(f"  sigma = {sigma*1e6:.1f} um: avg Rydberg = {mean_r:.2f} +/- {std_r:.2f}")

Position disorder has two effects:

1. It changes the blockade graph: If atoms are shifted closer together, new blockade edges appear. If shifted apart, some edges disappear. The ground state of the perturbed system may differ from the ideal case.
2. It breaks translational symmetry: Even if the blockade graph is unchanged, the interaction strengths V_ij are no longer uniform. This can split degeneracies in the ground state manifold and affect the dynamics of the adiabatic sweep.

For Aquila, typical position accuracy is well below 0.5 micrometers, so this effect is small for lattice spacings above 5 micrometers. But for tightly packed arrays or problems requiring precise interaction ratios, disorder can be the dominant source of error.

Benchmarking Against Classical Solvers

To understand whether the Rydberg simulator provides any advantage, compare its MIS solutions to those from a simple classical greedy algorithm.

import numpy as np
import random

def greedy_mis(adj, n):
    """Find an independent set using the greedy algorithm.
    
    Strategy: iteratively add the vertex with the fewest neighbors
    (among remaining candidates), then remove it and its neighbors.
    """
    remaining = set(range(n))
    independent_set = []

    while remaining:
        # Pick the vertex with the fewest remaining neighbors
        best_v = min(remaining, key=lambda v: len([u for u in adj[v] if u in remaining]))
        independent_set.append(best_v)
        # Remove this vertex and all its neighbors
        to_remove = {best_v} | {u for u in adj[best_v] if u in remaining}
        remaining -= to_remove

    return sorted(independent_set)

def random_greedy_mis(adj, n, n_trials=100):
    """Run greedy MIS with random vertex ordering multiple times."""
    best_size = 0
    best_set = []
    for _ in range(n_trials):
        remaining = set(range(n))
        indep_set = []
        vertices = list(remaining)
        random.shuffle(vertices)
        for v in vertices:
            if v not in remaining:
                continue
            indep_set.append(v)
            remaining -= {v} | set(adj[v])
        if len(indep_set) > best_size:
            best_size = len(indep_set)
            best_set = sorted(indep_set)
    return best_set

# Build the blockade graph for our 3x3 grid
lattice_spacing = 5.5e-6
positions = []
for row in range(3):
    for col in range(3):
        positions.append((col * lattice_spacing, row * lattice_spacing))

r_b = 8.5e-6  # blockade radius
adj = {i: [] for i in range(9)}
for i in range(9):
    for j in range(i + 1, 9):
        dx = positions[i][0] - positions[j][0]
        dy = positions[i][1] - positions[j][1]
        dist = np.sqrt(dx**2 + dy**2)
        if dist < r_b:
            adj[i].append(j)
            adj[j].append(i)

# Classical solutions
greedy_result = greedy_mis(adj, 9)
random_greedy_result = random_greedy_mis(adj, 9, n_trials=1000)

print("Classical greedy MIS:", greedy_result, f"(size {len(greedy_result)})")
print("Best random greedy MIS:", random_greedy_result, f"(size {len(random_greedy_result)})")

# Compare to AHS results
# Find the most common configuration from the AHS simulation
from collections import Counter

config_counts = Counter()
for shot in result.measurements:
    config = tuple(1 if s == 0 else 0 for s in shot.post_sequence)
    config_counts[config] += 1

top_configs = config_counts.most_common(5)
print("\nTop 5 AHS configurations (1=Rydberg, 0=ground):")
for config, count in top_configs:
    n_rydberg = sum(config)
    # Check if this is a valid independent set
    rydberg_sites = [i for i, c in enumerate(config) if c == 1]
    valid = True
    for i, a in enumerate(rydberg_sites):
        for b in rydberg_sites[i+1:]:
            if b in adj[a]:
                valid = False
                break
        if not valid:
            break
    status = "valid IS" if valid else "VIOLATION"
    print(f"  {config}: {count} shots, {n_rydberg} Rydberg atoms, {status}")

# For a 3x3 king graph, this is a trivially small problem.
# The real advantage of AHS appears for larger graphs (50+ vertices)
# where the MIS problem becomes genuinely hard for classical solvers.
print("\nNote: for a 3x3 grid, both classical and quantum methods find the")
print("optimal MIS easily. The potential quantum advantage appears for larger")
print("instances where the classical landscape has many local optima.")

For small graphs like the 3x3 grid, both classical and quantum methods find the optimal solution reliably. The hypothesis behind quantum optimization with Rydberg atoms is that for larger instances (50+ atoms), the quantum adiabatic dynamics can navigate the energy landscape more effectively than classical heuristics. This is an active area of research, and definitive quantum advantage for optimization has not yet been demonstrated.

From AHS Results to Optimization Solutions: Full Pipeline

Here is the complete pipeline for using Aquila to solve a graph optimization problem:

import numpy as np
from braket.ahs import (
    AtomArrangement, DrivingField, AnalogHamiltonianSimulation, TimeSeries
)
from braket.devices import LocalSimulator

# Step 1: Define the graph
# We use a unit disk graph: vertices are points in 2D, edges connect
# vertices within a threshold distance. This maps naturally to atom positions.

# Example: 9-vertex graph with specific connectivity
# Place atoms to create desired blockade edges
atom_positions_um = [
    (0.0, 0.0),    # vertex 0
    (5.5, 0.0),    # vertex 1: connected to 0
    (11.0, 0.0),   # vertex 2: connected to 1
    (0.0, 5.5),    # vertex 3: connected to 0
    (5.5, 5.5),    # vertex 4: connected to 1, 3 (and diagonals)
    (11.0, 5.5),   # vertex 5: connected to 2, 4
    (0.0, 11.0),   # vertex 6: connected to 3
    (5.5, 11.0),   # vertex 7: connected to 4, 6
    (11.0, 11.0),  # vertex 8: connected to 5, 7
]

# Step 2: Create the register (convert to meters)
register = AtomArrangement()
for x_um, y_um in atom_positions_um:
    register.add([x_um * 1e-6, y_um * 1e-6])

# Step 3: Build and run the AHS program
t_total = 4e-6
omega_drive = 2 * np.pi * 2.5e6
t_ramp = 0.5e-6

amp = TimeSeries()
amp.put(0.0, 0.0)
amp.put(t_ramp, omega_drive)
amp.put(t_total - t_ramp, omega_drive)
amp.put(t_total, 0.0)

det = TimeSeries()
det.put(0.0, -2 * np.pi * 5e6)
det.put(t_total, 2 * np.pi * 5e6)

ph = TimeSeries()
ph.put(0.0, 0.0)
ph.put(t_total, 0.0)

drive = DrivingField(amplitude=amp, detuning=det, phase=ph)
program = AnalogHamiltonianSimulation(register=register, hamiltonian=drive)

local_sim = LocalSimulator("braket_ahs")
result = local_sim.run(program, shots=1000).result()

# Step 4: Decode solutions
# Build the blockade graph for validation
r_b = 8.5e-6
adj = {i: [] for i in range(9)}
positions_m = [(x * 1e-6, y * 1e-6) for x, y in atom_positions_um]
for i in range(9):
    for j in range(i + 1, 9):
        dx = positions_m[i][0] - positions_m[j][0]
        dy = positions_m[i][1] - positions_m[j][1]
        if np.sqrt(dx**2 + dy**2) < r_b:
            adj[i].append(j)
            adj[j].append(i)

# Step 5: Extract and validate independent sets
valid_solutions = []
for shot in result.measurements:
    # Decode: post_sequence 0 = Rydberg (in the set), 1 = ground (not in the set)
    rydberg_sites = [i for i, s in enumerate(shot.post_sequence) if s == 0]

    # Validate: no two sites in the set should share a blockade edge
    is_valid = True
    for idx_a in range(len(rydberg_sites)):
        for idx_b in range(idx_a + 1, len(rydberg_sites)):
            a, b = rydberg_sites[idx_a], rydberg_sites[idx_b]
            if b in adj[a]:
                is_valid = False
                break
        if not is_valid:
            break

    if is_valid:
        valid_solutions.append(tuple(rydberg_sites))

# Find the best (largest) valid independent set
from collections import Counter
solution_counts = Counter(valid_solutions)
best_solution = max(solution_counts.keys(), key=lambda s: (len(s), solution_counts[s]))

print(f"Total shots: {len(result.measurements)}")
print(f"Valid independent sets: {len(valid_solutions)} ({len(valid_solutions)/len(result.measurements):.1%})")
print(f"Best solution: vertices {best_solution} (size {len(best_solution)})")
print(f"  Appeared {solution_counts[best_solution]} times")

# Print the solution on the grid
print("\nSolution on grid (R = in set, . = not in set):")
for row in range(3):
    row_str = ""
    for col in range(3):
        idx = row * 3 + col
        row_str += "R " if idx in best_solution else ". "
    print(f"  {row_str}")

Filtering for Atom Loading: Using pre_sequence

On real Aquila hardware, atom loading is probabilistic. Even though loading probability exceeds 99% per site, with 9 sites you expect roughly 9% of shots to have at least one missing atom. The pre_sequence field tells you which atoms were successfully loaded.

def filter_fully_loaded_shots(measurements):
    """Keep only shots where all atoms were successfully loaded."""
    full_shots = []
    partial_shots = 0
    for shot in measurements:
        if all(s == 1 for s in shot.pre_sequence):
            full_shots.append(shot)
        else:
            partial_shots += 1
    return full_shots, partial_shots

# On the local simulator, all atoms are always loaded.
# On real hardware, this filtering is essential.
full_shots, n_partial = filter_fully_loaded_shots(result.measurements)
print(f"Fully loaded shots: {len(full_shots)}")
print(f"Shots with missing atoms: {n_partial}")

# For more nuanced analysis, you can also handle partial loading
# by only analyzing the atoms that were present:
def analyze_with_partial_loading(measurements, register):
    """Analyze results while accounting for missing atoms."""
    positions = [(site[0], site[1]) for site in register]

    for shot in measurements[:5]:  # Show first 5 shots
        loaded = [i for i, s in enumerate(shot.pre_sequence) if s == 1]
        rydberg = [i for i, s in enumerate(shot.post_sequence) if s == 0]
        # Atoms that are both not loaded AND show 0 in post_sequence
        # are ambiguous: they could be lost (not Rydberg)
        ambiguous = [i for i in range(len(shot.pre_sequence))
                     if shot.pre_sequence[i] == 0 and shot.post_sequence[i] == 0]
        confident_rydberg = [i for i in rydberg if i not in ambiguous]

        print(f"  Loaded: {loaded}, Rydberg (confident): {confident_rydberg}, "
              f"Ambiguous: {ambiguous}")

analyze_with_partial_loading(result.measurements, register)

This distinction between Rydberg excitation and atom loss is important. When pre_sequence[i] == 0 (atom not loaded) and post_sequence[i] == 0, the zero in post_sequence indicates atom absence, not Rydberg excitation. Always cross-reference both arrays.

Running on Real Aquila Hardware

# To run on the actual QuEra Aquila device:
# from braket.aws import AwsDevice
#
# # Connect to Aquila
# aquila = AwsDevice("arn:aws:braket:us-east-1::device/qpu/quera/Aquila")
#
# # Check device availability
# print(f"Device status: {aquila.status}")
# print(f"Device properties: {aquila.properties}")
#
# # Submit the job (same AHS program as built above)
# task = aquila.run(ahs_program, shots=1000)
# print(f"Task ARN: {task.id}")
# print(f"Task status: {task.state()}")
#
# # Wait for results (this may take minutes to hours depending on queue)
# result = task.result()
#
# # Analysis code is identical to the local simulator case,
# # except you MUST filter for atom loading using pre_sequence.
# full_shots, n_partial = filter_fully_loaded_shots(result.measurements)
# print(f"Fully loaded: {len(full_shots)}, Partial: {n_partial}")

On real hardware, atom loading is probabilistic, so some sites may be empty. The pre_sequence field in each measurement indicates which sites were successfully loaded. Always filter your analysis to account for vacancies, as described in the previous section.

Common Mistakes and How to Avoid Them

Working with the Braket AHS API has several non-obvious pitfalls. Here are the most frequent mistakes and how to avoid them:

1. Misinterpreting post_sequence values

The encoding is counterintuitive: post_sequence[i] == 0 means the atom is in the Rydberg state (or was lost), and post_sequence[i] == 1 means the atom is in the ground state. This is the opposite of what most people expect. The convention comes from the detection scheme: the fluorescence imaging detects ground-state atoms (bright), while Rydberg atoms are dark.

# WRONG: treating 1 as Rydberg
n_rydberg_wrong = sum(1 for s in shot.post_sequence if s == 1)  # This counts GROUND state atoms

# CORRECT: treating 0 as Rydberg
n_rydberg_correct = sum(1 for s in shot.post_sequence if s == 0)  # This counts Rydberg atoms

2. Lattice spacing below 4 micrometers

The API rejects atom arrangements where any two atoms are closer than 4 micrometers. This constraint exists because the optical tweezers cannot reliably hold atoms closer than this.

# WRONG: spacing too small
# register.add([0, 0])
# register.add([3e-6, 0])  # 3 micrometers apart: will be rejected

# CORRECT: respect the minimum spacing
register.add([0, 0])
register.add([4e-6, 0])  # 4 micrometers: this is the minimum allowed

3. Ignoring pre_sequence on real hardware

On the local simulator, pre_sequence is always all ones (all atoms loaded). On real hardware, some atoms will be missing. If you skip the pre_sequence check, you will miscount Rydberg excitations: missing atoms show post_sequence == 0, which looks identical to Rydberg excitation.

# WRONG: analyzing post_sequence without checking pre_sequence
rydberg_sites = [i for i, s in enumerate(shot.post_sequence) if s == 0]

# CORRECT: only count sites that were loaded AND show Rydberg
rydberg_sites = [
    i for i in range(len(shot.post_sequence))
    if shot.pre_sequence[i] == 1 and shot.post_sequence[i] == 0
]

4. Running too few shots

With 9 atoms and 2^9 = 512 possible configurations, you need enough shots to resolve the distribution. A minimum of 1000 shots is recommended for 9-atom systems. For larger systems, the configuration space grows exponentially, but the physically relevant configurations (those respecting the blockade) are a much smaller subset.

5. Trusting the local simulator as a hardware model

The LocalSimulator("braket_ahs") solves the Schrodinger equation for the Rydberg Hamiltonian exactly. It does not model:

• Atom loading failures (pre_sequence is always all ones)
• Detection errors (post_sequence is sampled from the ideal quantum state)
• Laser phase noise
• Finite temperature effects
• Atom loss during evolution

For these reasons, local simulator results are an upper bound on the solution quality you can expect from real hardware. The simulator is useful for validating your pulse design and understanding the ideal physics, but do not treat its output as a prediction of hardware performance.

6. Omega not starting and ending at zero

The hardware requires that the Rabi frequency starts at zero and ends at zero. If you try to create a TimeSeries where amplitude.put(0.0, omega_max), the API will reject it. Always include ramp-up and ramp-down segments.

Key Takeaways

1. Analog Hamiltonian simulation on Aquila is programmed through atom positions and time-dependent laser pulses, not gates. The atom positions encode the problem structure through the distance-dependent blockade interaction, and the pulse shapes control the adiabatic evolution.

2. The Rydberg blockade radius determines the interaction graph. For Omega = 2*pi * 2.5 MHz, r_b ≈ 8.5 micrometers. Atoms closer than r_b cannot both be in the Rydberg state, creating an effective constraint graph. The relationship between lattice spacing and blockade radius is the primary design knob.

3. The adiabatic sweep maps Rydberg physics to combinatorial optimization. Sweeping Delta from negative to positive drives the system from the trivial all-ground state to the Maximum Independent Set of the blockade graph. The sweep rate determines the solution quality.

4. Local detuning enables weighted optimization problems. By assigning different detuning values to different atoms, you can encode vertex weights and solve Maximum Weight Independent Set.

5. Spatial correlations reveal the structure of the prepared state. The connected correlation function C(r) distinguishes random excitations from true blockade-ordered phases. Anti-correlation at short distances is the blockade signature.

6. Position disorder, atom loading, and decoherence are the main sources of error on real hardware. The local simulator does not model these effects, so always validate on hardware before drawing conclusions about quantum advantage.

7. Always check pre_sequence on real hardware to distinguish between Rydberg excitation and atom loss. Filter or correct for missing atoms before analyzing post_sequence.

8. The potential for quantum advantage in optimization is an open research question. For small systems, classical solvers easily match or exceed the quantum results. The hypothesis is that larger, more complex instances will eventually reveal a quantum advantage. Current research focuses on identifying the problem classes and parameter regimes where this advantage is most likely to emerge.