Skip to main content

View on GitHub

Open this notebook in GitHub to run it yourself
Quantum maps, the quantum analogs of classical maps, are simple yet powerful models of quantum dynamics. A typical quantum map describes a system undergoing free evolution that is periodically “kicked” by a position-dependent force. The resulting dynamics are discrete in time and are naturally formulated as a so-called Floquet system [1]. Most quantum maps are obtained by quantizing a known classical map, such as the tent map, the sawtooth map, or the kicked rotor, allowing for direct comparison between classical and quantum behavior. Over the years, these toy models have played an important role in understanding fundamental phenomena such as the correspondence principle, decoherence, entanglement spreading, and quantum dynamical chaos. In the context of quantum computing, these well-studied systems provide a convenient and nontrivial testbed for quantum algorithms and hardware. Their dynamics are rich enough to generate strong entanglement and chaotic behavior, yet structured enough to admit efficient implementations on quantum circuits. As a result, quantum maps have been used as controlled settings to investigate the effects of noise, imperfections, and decoherence in quantum processors, as well as to assess hardware performance. For example, quantum chaotic maps have been employed to study the impact of realistic noise models on quantum dynamics [2] and to benchmark quantum hardware through their sensitivity to errors and entanglement generation [3]. Despite their simple definition, quantum maps can exhibit chaotic dynamics with rapid entanglement growth, which makes classical simulation inefficient. In particular, it has been shown that these systems can be simulated exponentially faster on a quantum computer [4]. Since these systems are already well studied, this exponential speedup is not expected to translate into a scientific quantum advantage today. In practice, classical numerical simulations at sizes on the order of ~20 qubits are typically sufficient to capture and understand most of the phenomena currently of interest. In this notebook we focus on the quantum sawtooth map and the phenomenon of dynamical localization, whose validation on quantum hardware has been studied in several previous works (see Refs. in [1]). Below, we define the model and discuss its theoretical properties, as well as present execution results on an IonQ machine. This notebook illustrates how Classiq Qmod’s language and execution framework can be used to easily define, run, and observe quantum phenomena on a real quantum device.

​
Quantum Sawtooth Map and its Evolution

​

​
The system Hamiltonian

We start with the classical sawtooth map, written in terms of the action-angle variables (p,q)(p,q) Hcl(t)=p22L−Kq22∑iδ(t−iτ),q mod 2π,\large H_{\rm cl}(t) = \frac{p^2}{2L}-K\frac{q^2}{2}\sum_{i}\delta(t-i\tau), \quad q \bmod 2\pi, where LL is the moment of inertia, and KK and τ\tau are the kicking strength and kicking period, respectively. We now skip directly to the quantum model; the full derivation is provided in the Technical Notes at the end of this notebook. The derivation includes a straightforward quantization procedure, together with several redefinitions and normalizations. For the quantum map, we consider a finite Hilbert space of dimension NN, with momentum and position operators p^\hat{p} and q^\hat{q} defined by p^∣p⟩=p∣p⟩,q^∣q⟩=q∣q⟩,p,q∈{N2,…,−1,0,1,…,N−12}.\hat{p}|p\rangle = p|p\rangle,\qquad \hat{q}|q\rangle = q|q\rangle, \qquad p,q \in \left\{\frac{N}{2}, \dots,-1,0,1, \dots,\frac{N-1}{2}\right\}. The Hamiltonian of the quantum sawtooth map reads: HN(t)ℏ=2π(p^22N−Kq^22N∑iδ(t−i)),\large \frac{H_N(t)}{\hbar} = 2\pi \left(\frac{\hat{p}^2}{2N}-K\frac{\hat{q}^2}{2N}\sum_{i}\delta(t-i)\right), where ℏ\hbar is a dimensionless Planck constant (which can be related to the physical one by a factor of L/τL/\tau). As part of the quantization procedure, the product ℏN\hbar N is held constant. Consequently, the parameter controlling quantum effects is the Hilbert space dimension NN, with quantum fluctuations scaling as ℏ∼1/N\hbar \sim 1/N.

​

​
The time evolution

A key practical advantage of quantum maps is the simplicity of their implementation on a quantum computer. Because the kicks are modeled as delta functions, the time-evolution operator factorizes into a sequence of evolutions generated by operators that are diagonal in complementary bases. In particular, for the sawtooth map considered in this notebook, the implementation does not require approximation overhead from product formulas or Chebyshev expansions, making it especially well suited for near-term quantum devices. The unitary evolution of the system between kick tt and t+1t+1, known as the Floquet operator, is given by: ∣ψt+1⟩=Ukick∣ψt⟩,Ukick≡exp⁡[−i2πNp^22]exp⁡[iK2πNq^22].\large |\psi_{t+1}\rangle = U_{\rm kick}|\psi_{t}\rangle,\qquad U_{\rm kick} \equiv \exp\left[-i\frac{2\pi}{N} \frac{\hat{p}^2}{2}\right] \exp\left[iK\frac{2\pi}{N} \frac{\hat{q}^2}{2}\right]. (This expression is exact and follows from modeling the kicks as delta-function impulses.) The evolution of the system from an initial state after TT kicks is therefore UT−kicks=UkickT=(exp⁡[−i2πNp^22]exp⁡[iK2πNq^22])T.\large U_{T - {\rm kicks}} = U^T_{\rm kick} = \left(\exp\left[-i\frac{2\pi}{N} \frac{\hat{p}^2}{2}\right] \exp\left[iK\frac{2\pi}{N} \frac{\hat{q}^2}{2}\right]\right)^T. The kinetic term is diagonal in the momentum basis, while the kicking term is diagonal in the position basis. Since the two bases are related by a (discrete) Fourier transform, the Floquet operator can be efficiently simulated by alternating between these representations.

​

​
Implementation in Qmod

We work with a Hilbert space of nn qubits, corresponding to a dimension N=2nN=2^n. The time evolution of the model can be written in Qmod using only a few lines of code, taking advantage of the phase assignment construct. 
from classiq import *
from classiq.qmod.symbolic import pi


@qfunc
def p_space_evolution(p: QNum):
    phase(-(pi / 2**p.size) * p**2)


@qfunc
def q_space_evolution(kick_mag: CReal, p: QNum):
    within_apply(lambda: qft(p), lambda: phase((kick_mag * pi / (2**p.size)) * p**2))


@qfunc
def single_step(kick_mag: CReal, p: QNum):
    q_space_evolution(kick_mag, p)
    p_space_evolution(p)

​
Dynamical Localization

Depending on the non-dimensional parameters NN and KK, the quantum system exhibits different dynamical regimes and a variety of phenomena, such as dynamical localization, anomalous diffusion, and a positive quantum Lyapunov exponent. In this notebook we focus on the regime 0<K<10<K<1, in which dynamical localization can be observed (we refer to the literature, e.g., Ref. [1], for the full plethora of possible behaviors). For 0<K<10<K<1,starting from a single excited momentum eigenstate ∣p(0)⟩=∣p0⟩|p(0)\rangle = |p_0\rangle, the wavefunction initially diffuses classically in momentum space, until quantum effects become relevant. This occurs at the Heisenberg time, defined as the time required to resolve the energy levels of the system. Beyond this time scale, quantum interference suppresses diffusion and leads to localization of the wavefunction over long times, ∣⟨p∣ψ⟩∣2∼1ℓe−2∣p−p0∣/ℓ,\left|\langle p|\psi\rangle\right|^2 \sim \frac{1}{\ell}e^{-2|p-p_0|/\ell}, with localization length ℓ≈3.3K5/2N2/(2π)2\ell\approx 3.3 K^{5/2}N^2/(2\pi)^2. Since we consider a finite Hilbert space, localization persists only for a finite time and is observable provided that the localization length is smaller than the system size, ℓ<N\ell<N. This condition is satisfied as long as K<Kloc≈πN−2/5.K<K_{\rm loc} \approx \pi N^{-2/5}.

​
Example: Localization in a 3 qubits System

We choose some hyper parameters, looking at a system on 3 qubits. 
import numpy as np

NUM_QUBITS = 3
hilbert_size = 2**NUM_QUBITS
k_loc = np.pi * hilbert_size ** (-2 / 5)
print(f"Kicking magnitude must be smaller than K={k_loc}")

KICK_MAG = 0.1
print(f"Taking K = {KICK_MAG}")
loc_len = 3.3 * KICK_MAG ** (5 / 2) * (hilbert_size / 2 * np.pi) ** (2)
assert (
    loc_len < hilbert_size
), "Localization length {loc_len} must be smaller than system size"
print(f"Localization length:  {loc_len}")
Output:

Kicking magnitude must be smaller than K=1.3674576271147798
  Taking K = 0.1
  Localization length:  1.64791067824756

​

​
Model definition

Defining a model with a parametric number of total time (number of kicks TT): 
@qfunc
def main(num_kicks: CInt, kick_mag: CReal, p: Output[QNum[NUM_QUBITS, SIGNED, 0]]):
    allocate(p)
    p ^= -2
    power(num_kicks, lambda: single_step(KICK_MAG, p))


qprog = synthesize(main)
show(qprog)
Output:

Quantum program link: https://platform.classiq.io/circuit/39hMn1YM3alzTTAUqzeRgi7MHR4
Screenshot 2026-02-05 at 22.24.07.png 
depth = qprog.transpiled_circuit.depth
cx_counts = qprog.transpiled_circuit.count_ops["cx"]
print(f"depth for one iteration {depth}")
print(f"cx counts for one iteration {cx_counts}")
Output:
depth for one iteration 48
  cx counts for one iteration 30

​

​
Execution

We execute the quantum program for different total numbers of kicks in order to examine the time evolution of the initial wavefunction. 
NUM_TIMES = 4
timesteps = [2 * i for i in range(1, NUM_TIMES + 1)]


print(f"======== Running for {len(range(1,NUM_TIMES+1))} different times =========")
print(f"times : {timesteps}")
print(f"depths per power execution {[depth*2*i for i in range(1,NUM_TIMES+1)]}")
cxs = [cx_counts * 2 * i for i in range(1, NUM_TIMES + 1)]
print(f"cx-counts per power execution  {cxs}")
print(f"Total cx counts {sum(cxs)}")
Output:
======== Running for 4 different times =========
  times : [2, 4, 6, 8]
  depths per power execution [96, 192, 288, 384]
  cx-counts per power execution  [60, 120, 180, 240]
  Total cx counts 600
First, we run the circuit on a simulator to obtain a reference for the expected dynamics. 
NUM_SHOTS = 1000


backend_preferences = ClassiqBackendPreferences(
    backend_name="simulator",
)
execution_prefs = ExecutionPreferences(
    backend_preferences=backend_preferences, num_shots=NUM_SHOTS
)

with ExecutionSession(qprog, execution_prefs) as es:
    results_simulator = es.batch_sample([{"num_kicks": t} for t in timesteps])
When running on hardware, we submit a job and store its ID. The results can then be retrieved once the execution is complete. Here is an example of how to run on IonQ machines, using run_through_classiq. 
import datetime
import os

BACKEND_NAME = "qpu.forte-1"
RUN_ON_HARDWARE = False

if RUN_ON_HARDWARE:
    # Define backend preferences

    backend_preferences = IonqBackendPreferences(
        backend_name=BACKEND_NAME, run_through_classiq=True, error_mitigation=False
    )

    prefix = f"{BACKEND_NAME}_{NUM_SHOTS}"
    execution_prefs = ExecutionPreferences(
        backend_preferences=backend_preferences, num_shots=NUM_SHOTS
    )
    # Submit a job
    with ExecutionSession(qprog, execution_prefs) as es:
        job = es.submit_batch_sample([{"num_kicks": t} for t in timesteps])
        job_ID = job.id
        print(f"Job ID: {job_ID}")
        timestamp = datetime.datetime.now().strftime("%Y%m%d_%H%M%S")
        filename = f"{prefix}_jobID_{timestamp}.txt"

        with open(filename, "w") as f:
            f.write(job_ID)


# Check job status
if RUN_ON_HARDWARE:
    with open(filename, "r") as f:
        job_ID = f.read()

    restored_job = ExecutionJob.from_id(job_ID)
    status = restored_job.status
    print(f"Job status: {status}")

    if status == "COMPLETED":
        res = restored_job.get_batch_sample_result()
        dirname = f"data_{prefix}_jobID_{timestamp}"
        os.makedirs(dirname, exist_ok=True)

        for i, df in enumerate(res):
            df = res[i].dataframe
            df.to_csv(f"{dirname}/df_{timesteps[i]}kicks.csv", index=False)
Below, we load pre-run results obtained on IonQ and IBM machines, and compare them to the simulator results: 
dirname_ionq = "data_qpu.forte-1_1000_jobID_20260209_150802"


def sort_and_fill_df(df):
    half = 2 ** (NUM_QUBITS - 1)
    p_full = np.arange(-half, half + 1)

    df_sorted = (
        df.set_index("p")  # make p the index
        .reindex(p_full)  # add missing p values
        .fillna({"probability": 0.0})  # fill missing probability with 0
        .reset_index()
        .rename(columns={"index": "p"})
        .sort_values("p")
    )

    return df_sorted


import matplotlib.pyplot as plt
import pandas as pd

fig, axes = plt.subplots()

# Simulator
for i in range(len(timesteps)):
    df = results_simulator[i].dataframe
    df_sorted = sort_and_fill_df(df)
    cl_plot = axes.plot(
        df_sorted["p"], df_sorted["probability"], ":o", color="gray", linewidth=4
    )

axes.tick_params(axis="both", labelsize=14)
pl_objs = []
for i in range(len(timesteps)):
    df = pd.read_csv(f"{dirname_ionq}/df_{timesteps[i]}kicks.csv")
    df_sorted = sort_and_fill_df(df)
    pl_objs += axes.plot(df_sorted["p"], df_sorted["probability"], "-o", linewidth=2)

axes.legend(
    cl_plot + pl_objs,
    [f"simulator, T= {', '.join(map(str, timesteps))}"]
    + [f"qpu, T = {timesteps[i]}" for i in range(len(timesteps))],
    fontsize=14,
)
get_backend_name = dirname_ionq.split("data_", 1)[1].split(f"_{NUM_SHOTS}", 1)[0]
axes.set_title(f"{get_backend_name}", fontsize=16)
axes.set_xlabel(r"$p$", fontsize=14)
axes.set_ylabel(r"$P(p)$", fontsize=14)

plt.tight_layout()
plt.show()
output 
for i in range(len(timesteps)):
    df = results_simulator[i].dataframe
    assert float(df.loc[df["p"] == -2, "probability"].iloc[0]) > 0.85

​
Technical Notes

Below we derive the time-evolution of the quantum sawtooth map. We start with the classical map and go through several renormalizations and redefinitions. The classical system is written for the so-called action-angle variables (p,q)(p,q) Hcl(t)=p22L−Kq22∑iδ(t−iτ),q mod 2π,\large H_{\rm cl}(t) = \frac{p^2}{2L}-K\frac{q^2}{2}\sum_{i}\delta(t-i\tau), \quad q \bmod 2\pi, where LL is the momentum of inertia, and KK and τ\tau are the kicking strength and period, respectivaly. We move to a non-dimentional problem, by taking t→tτt\rightarrow t\tau, and p→pL/τp\rightarrow pL/\tau Hcl(t)=p22−Kq22∑iδ(t−i),q mod 2π,\large H_{\rm cl}(t) = \frac{p^2}{2}-K\frac{q^2}{2}\sum_{i}\delta(t-i), \quad q \bmod 2\pi, where we also define K→Kτ/LK\rightarrow K\tau/L. Next, we quantize the system, defining the operators p^\hat{p} and q^\hat{q} on a finite Hilbert space of size NN with p^∣n⟩=ℏ∣n⟩,[p^,q^]=iℏ,\hat{p}|n\rangle = \hbar|n\rangle, \quad [\hat{p},\hat{q}]=i\hbar, p=N2,…,−1,0,1,…,N−12,p = \frac{N}{2}, \dots,-1,0,1, \dots,\frac{N-1}{2}, q/2π=−12,⋯−1N,0,1N…N−12N.q/2\pi= -\frac{1}{2},\dots -\frac{1}{N},0,\frac{1}{N}\dots \frac{N-1}{2N}. where ℏ\hbar is the dimensionless Planck’s constant (can be related the physical one by a factor of L/τL/\tau). The quantized Hamiltonian reads H(t)=p^22−Kq^22∑iδ(t−i).\large H(t) = \frac{\hat{p}^2}{2}-K\frac{\hat{q}^2}{2}\sum_{i}\delta(t-i). ` As a final step, we consider a fixed phase-space area, defined according to ℏN=(2π)⋅Constant=2π,\hbar N = (2\pi)\cdot \text{Constant}=2\pi, (for the semi-classical limit this is just the Bohr-Sommerfeld approximation), and redefine: p^→ℏp,q→2πNq^.\hat{p}\rightarrow \hbar{p}, q\rightarrow \frac{2\pi}{N}\hat{q}. This final transformation results in the Hamiltonian HN(t)ℏ=2π(p^22N−Kq^22N∑iδ(t−i)).\large \frac{H_N(t)}{\hbar} = 2\pi \left(\frac{\hat{p}^2}{2N}-K\frac{\hat{q}^2}{2N}\sum_{i}\delta(t-i)\right). `

​
References

[1]: M. S. Rudner, N. H. Lindner. “The Floquet Engineer’s Handbook”. arXiv:2003.08252 [cond-mat.mes-hall] (2020) [2]: M. D. Porter, I. Joseph, “Impact of dynamics, entanglement, and Markovian noise on the fidelity of few-qubit digital quantum simulation”. Journal of Plasma Physics 91.1 (2025) [3]: Andrea Pizzamiglio, Su Yeon Chang, Maria Bondani, Simone Montangero, Dario Gerace, Giuliano Benenti. “Dynamical localization simulated on actual quantum hardware.” Entropy 23, 654 (2021). [4]: Georgeot, B. and Shepelyansky, D.L., 2001. “Exponential gain in quantum computing of quantum chaos and localization”. Physical Review Letters, 86(13), p.2890.