Sequences and Simulations

Sequence

Container for a time-ordered list of low-level control instructions with a fixed time step dt.

A Sequence behaves like a vector of AbstractInstruction objects: it supports indexing, iteration, length, and push!, and is the primary container passed to play when simulating instruction-level dynamics.

Fields

• instructions::Vector{AbstractInstruction} – ordered list of instructions in the sequence
• dt – base time step used to interpret instruction durations (typically a Float64 in seconds or an application-specific time unit)

Sequences are usually constructed via the convenience constructor

Sequence(dt::Float64)

which creates an empty sequence with the given time step.

Sequence(dt::Float64)

Create an empty Sequence with time step dt and no instructions.

This is the recommended constructor to start building a new instruction sequence, e.g.

seq = Sequence(1.0e-9) # 1 ns base time step

Base.eltype(::Type{Sequence})

Return the element type of a Sequence, which is AbstractInstruction.

Base.getindex(seq::Sequence, i...)

Index into the underlying instruction list of seq.

This allows Sequence to be used like a vector of AbstractInstruction, e.g. seq[1] returns the first instruction in the sequence.

Base.iterate(seq::Sequence, state...)

Iterate over the instructions in seq.

Enables use of for inst in seq and other iterator-based patterns, treating Sequence as a collection of AbstractInstruction objects.

Base.length(seq::Sequence)

Return the number of instructions stored in seq.

Base.push!(seq::Sequence, inst::AbstractInstruction)

Append an instruction inst to the end of seq.

Returns the modified Sequence, allowing idioms such as

push!(seq, Wait(10e-9))
push!(seq, Pulse(global_coupling, t))

@sequence seq begin 
    ...
end

Macro for building an instruction Sequence from imperative-style code while preserving native Julia control flow.

@sequence walks the body of the begin ... end block and rewrites each function call into a push!(seq, ...). Loops, conditionals, and other control-flow constructs are left intact and execute as normal, so you can generate instruction patterns programmatically:

@sequence seq begin
    Wait(1µs)
    for i in 1:3
        Pulse(coupling, t)
        Wait(i * 500ns)
    end
    if true
        Wait(10ns)
    end
end

AbstractInstruction

Abstract supertype for all low-level control instructions (pulses, switching, waits, tweezer moves, ramps, etc.).

Pulse(couplings::Vector{<:Switchable}, duration, ampl, amplitudes)

Apply laser coupling for a fixed duration or with a shaped amplitude profile.

Fields

• couplings::Vector{Switchable} – list of Switchable objects (e.g. GlobalCoupling{NLevelAtom}) to be modulated.
• duration::Float64 – pulse duration in seconds. The compiled pulse will span exactly this duration.
• ampl::ComplexF64 – overall coupling strength (relative to the default). Can be complex. It scales the entire amplitude shape.
• amplitudes::Vector{ComplexF64} – vector of complex amplitudes defining the pulse shape. If empty, the pulse has constant amplitude ampl over duration. If non‑empty, the shape defined by amplitudes is resampled and interpolated over the full duration duration (using third‑order Lagrange interpolation).

On(couplings)

Turn on one or more couplings.

• couplings is a vector of Switchable objects to be switched on.

Off(couplings)

Turn off one or more couplings.

• couplings is a vector of Switchable objects to be switched off.

Wait(duration)

Idle period with no explicit control actions.

• duration is the wait time (seconds).

MoveCol(tweezers, cols, delta, duration; sweep = :min_jerk)

Move one or more columns of tweezers simultaneously.

• cols can be a single integer or a vector of integers specifying the columns.
• delta is the frequency shift (Hz).
• duration is the total time for the move (seconds).
• sweep is the sweep profile: built-in symbol (:linear, :cosine, :min_jerk) or custom function s -> f(s) mapping [0,1] → [0,1].

MoveRow(tweezers, rows, delta, duration; sweep = :min_jerk)

Move one or more rows of tweezers simultaneously.

• rows can be a single integer or a vector of integers specifying the rows.
• delta is the frequency shift (Hz).
• duration is the total time for the move (seconds).
• sweep is the sweep profile: built-in symbol (:linear, :cosine, :min_jerk) or custom function s -> f(s) mapping [0,1] → [0,1].

RampCol(tweezers, cols, final_amplitude, ramp_time)

Ramp the amplitude of one or more columns of tweezers simultaneously.

• cols can be a single integer or a vector of integers specifying the columns.
• final_amplitude is the value to ramp each column to (scalar or vector).
• ramp_time is the total duration of the ramp (seconds).

AmplCol(tweezers, col, ampl)

Set the relative amplitude of a single tweezer column.

• col is the column index.
• ampl is the new relative amplitude.

AmplRow(tweezers, row, ampl)

Set the relative amplitude of a single tweezer row.

• row is the row index.
• ampl is the new relative amplitude.

FreqCol(tweezers, col, freq)

Set the frequency of a single tweezer column.

• col is the column index.
• freq is the absolute frequency (Hz).

FreqRow(tweezers, row, freq)

Set the frequency of a single tweezer row.

• row is the row index.
• freq is the absolute frequency (Hz).

play(sys::System, seq::Sequence; 
        initial_state=sys.initial_state, 
        rng=Random.default_rng(), kwargs...) -> NamedTuple

Execute a quantum simulation by compiling and running a pulse sequence on a system.

This is the high-level entry point for running simulations. It automatically handles compilation, state initialization, and execution in a single call. For performance-critical workflows with repeated executions, consider using compile followed by play(::SimulationJob, ::System) to avoid recompilation overhead.

Arguments

• sys::System: System specification containing atoms, beams, operators, and detector configurations
• seq::Sequence: Time-ordered instruction sequence (pulses, moves, ramps, waits) with timestep dt
• initial_state: Initial quantum state specification (required for quantum systems). Can be:
  • AbstractVector: Basis-ordered state vector or density matrix
  • Tuple: Collection of basis levels, e.g., (g, g, e) for three atoms
  • AbstractLevel: Single level for uniform initialization
  • Default: sys.initial_state

Keyword Arguments

• shots::Int = 1: Number of Monte Carlo trajectory shots to execute
• density_matrix::Bool = false: Use density matrix formalism if true, statevector if false
• savefinalstate::Bool = false: Include final quantum states in output (increases memory usage)
• rng::AbstractRNG = Random.default_rng(): Random number generator for reproducible simulations
• Additional kwargs are treated as parameter values for resolving Deferred objects and other parametric components in the system and sequence

Returns

Returns a NamedTuple with the following fields:

• detectors::Dict{String, Array}: Detector measurement outputs. Format depends on shots:
  • Single shot (shots=1): Dict{String, Vector} for 1D detectors, Dict{String, Matrix} for multi-dimensional
  • Multiple shots (shots>1): [n_times × shots] for 1D detectors, [n_times × n_dims × shots] for multi-dimensional
• times::Vector{Float64}: Global time points at which measurements were recorded (starts at dt, not zero)
• final_states::Vector: Final quantum state after evolution for each shot (only if savefinalstate=true)

Notes

• For quantum systems, initial_state must be specified in sys.initial_state or overridden via this argument
• Each shot reinitializes atomic velocities/positions with fresh randomness
• Detector measurements occur at the end of each timestep, not the beginning
• Multi-shot simulations use parallel execution when shots ≥ 4 and Threads.nthreads() > 1
• Classical systems (no quantum state) skip quantum evolution and only simulate atomic motion

Examples

Single-shot quantum simulation

```julia using AtomTwin

Define system with ground and excited states

g, e = AtomTwin.Level(:g), AtomTwin.Level(:e) atoms = [Atom(position=[0.0, 0.0, 0.0], levels=[g, e])] basis = Basis([g, e])

Create Rabi pulse and detector

Ω = RabiField(amplitude=2π*1.0, detuning=0.0) detector = PopulationDetector("P_g", level=g)

sys = System(atoms, [], [g], basis, [Ω], [detector]) seq = Sequence(dt=0.01) push!(seq, Pulse(Ω, duration=1.0))

Run simulation

result = play(sys, seq; initialstate=g, shots=1) result.detectors["Pg"] # Vector of ground state populations vs time

Multi-shot Monte Carlo simulation

Same system as above, run 100 trajectories with quantum jumps

jump = Jump(rate=0.1, source=e, target=g) sys_jumps = System(atoms, [], [g], basis, [Ω, jump], [detector])

result = play(sysjumps, seq; initialstate=g, shots=100, savefinalstate=true) result.detectors["Pg"] # Matrix: [ntimes × 100] meanpopulation = mean(result.detectors["Pg"], dims=2) # Average over shots

Parametric simulation with deferred values

Define amplitude as a parameter

Ωparam = Deferred(:amplitude) pulse = Pulse(RabiField(amplitude=Ωparam, detuning=0.0), duration=1.0)

Run with specific parameter value

result = play(sys, seq; initial_state=g, amplitude=2π*2.0)

compile(atoms, inst, dt; resolve_target = identity) -> (modifiers, tspan)

Lower an inert instruction spec inst into concrete modifiers and time spans for a single instruction.

Common contract:

• Inputs:

  • atoms: current atom state (vector or model-specific container)
  • inst: instruction spec (e.g. MoveCol, Wait, RampCol, Pulse, …)
  • dt: time step (seconds, Float64)
  • resolve_target: function that binds objects captured by specs to runtime instances

• Output:

  • (modifiers, tspan) where modifiers is a vector of per-segment modifiers and tspan is the segment time vector

• Time convention:

  • tspan is local to the segment, typically collect(dt:dt:duration)
  • Downstream code accumulates absolute time across segments

• Notes:

  • Off, Pulse instructions set the amplitude of the target coupling to 0.0 at the end of the instruction

compile(atoms, inst::MoveRow, dt; resolve_target = identity)

Lower a MoveRow instruction into position modifiers that move tweezers in the specified rows along y according to the chosen sweep profile.

compile(atoms, inst::MoveCol, dt; resolve_target = identity)

Lower a MoveCol instruction into position modifiers that move tweezers in the specified columns along x according to the chosen sweep profile.

compile(atoms, inst::Wait, dt; resolve_target = identity)

Lower a Wait instruction into an idle time segment. No modifiers are produced.

compile(atoms, inst::RampRow, dt; resolve_target = identity)

Lower a RampRow instruction into amplitude modifiers that linearly ramp the amplitudes of tweezers in the selected columns over inst.ramp_time.

compile(atoms, inst::RampCol, dt; resolve_target = identity)

Lower a RampCol instruction into amplitude modifiers that linearly ramp the amplitudes of tweezers in the selected columns over inst.ramp_time.

compile(atoms, inst::AmplCol, dt; resolve_target = identity)

Lower an AmplCol instruction into single-step amplitude modifiers for all tweezers in the given column, enforcing row–column factorization.

compile(atoms, inst::AmplRow, dt; resolve_target = identity)

Lower an AmplRow instruction into single-step amplitude modifiers for all tweezers in the given row, enforcing row–column factorization.

compile(atoms, inst::FreqCol, dt; resolve_target = identity)

Lower a FreqCol instruction into position modifiers that update the effective position (e.g. optical frequency) of all beams in the given column.

compile(atoms, inst::FreqRow, dt; resolve_target = identity)

Lower a FreqRow instruction into position modifiers that update the effective position (e.g. optical frequency) of all beams in the given row.

compile(atoms, inst::On, dt; resolve_target = identity)

Lower an On instruction into a single-step amplitude modifier that sets the amplitude of the given couplings to 1 at t = 0.

compile(atoms, inst::Off, dt; resolve_target = identity)

Lower an Off instruction into a single-step amplitude modifier that sets the amplitude of the given couplings to 0 at t = 0.

compile(atoms, inst::Parallel, dt; resolve_target = identity)

Compile several instructions to execute in parallel over a common time interval equal to the longest internal instruction. Updates that exceed length(m.vals) for modifier m are effectively ignored by the modifier.

compile(system::System, sequence::Sequence; initial_state=nothing, density_matrix=false) -> SimulationJob

Compile a System and Sequence into an executable SimulationJob with preallocated single-shot storage.

Arguments

• system: System specification (atoms, beams, nodes, detectors)
• sequence: Pulse sequence to execute (instruction list with timestep dt)
• initial_state: Initial quantum state (required for quantum systems)
• density_matrix: Use density matrix formalism if true (default: false)
• Additional keyword arguments are treated as parameter overrides (e.g. Ω = 2π*1e6)

Returns

• SimulationJob ready for execution, containing single-shot detector output buffers

Notes

• Compiles all DAG nodes: samples parameter values and updates fields in-place
• Detectors are automatically type-specialized to avoid dynamic dispatch
• Detector outputs are preallocated views into storage, avoiding allocations during simulation
• Each call to compile() creates a single-shot job. Multi-shot execution in play() uses thread-local copies of this job, with results aggregated into output matrices.

recompile!(job::SimulationJob, sys::System; kwargs...)

Reinitialize a SimulationJob for a new Monte Carlo trajectory.

Updates all DAG node outputs in-place (re-sampling parameter values and noise), reinitializes atom velocities, resets the quantum state, and zeroes detector outputs.

Additional keyword arguments override parameter values (same as compile).

Thread Safety

Safe to call on thread-local job copies (deepcopy(job)). MUST NOT be called on shared job objects. The sys argument may be shared across threads.

evolve!(atoms, system, tspan; frozen = false, kwargs...)

High-level entry point for pure classical motion of atoms atoms in a composite system containing beams and possibly other components.

All atoms are initialized in their ground level population, and Newtonian equations of motion are integrated by newton.

evolve!((rho, atoms), tspan; fields, jumps, beams, frozen, kwargs...)

Evolve a density matrix rho and atomic positions atoms over tspan using either a frozen-atom quantum master equation or a semiclassical master equation with classical trajectories.

evolve!((psi, atoms), tspan; fields, jumps, beams, frozen, kwargs...)

Evolve a state vector psi and atomic positions atoms over tspan using either deterministic Schrödinger dynamics or a wavefunction Monte Carlo trajectory, with optional semiclassical motion.