Library Reference

AbstractDrive

Abstract supertype for Hamiltonian drive terms.

A drive pairs a Hermitian matrix H with a scalar coefficient that depends on the control vector u. The full Hamiltonian is:

H(u, t) = H_drift + Σ_d drive_coeff(d, u) * d.H

Subtypes must implement:

• drive_coeff(d, u) — compute the scalar coefficient
• drive_coeff_jac(d, u, j) — compute ∂coeff/∂u_j
• active_controls(d) — return indices where ∂coeff/∂u_j can be nonzero

See LinearDrive and NonlinearDrive.

AbstractQuantumSystem

Abstract type for defining systems.

CompositeQuantumSystem <: AbstractQuantumSystem

A composite quantum system consisting of multiple subsystems with optional coupling terms.

Composite systems represent multiple quantum subsystems (e.g., multiple qubits or oscillators) that may be coupled together. Each subsystem's Hamiltonians are automatically lifted to the full tensor product space, and subsystem drives are appended to any coupling drives.

Fields

• H::Function: The total Hamiltonian function: (u, t) -> H(u, t)
• G::Function: The isomorphic generator function: (u, t) -> G(u, t)
• H_drift::SparseMatrixCSC{ComplexF64, Int}: The total drift Hamiltonian including subsystem drifts and couplings
• H_drives::Vector{SparseMatrixCSC{ComplexF64, Int}}: All drive Hamiltonians (coupling drives + subsystem drives)
• drive_bounds::Vector{Tuple{Float64, Float64}}: Drive amplitude bounds for each control
• n_drives::Int: Total number of control drives
• levels::Int: Total dimension of the composite system (product of subsystem dimensions)
• subsystem_levels::Vector{Int}: Dimensions of each subsystem
• subsystems::Vector{QuantumSystem}: The individual quantum subsystems

See also Lifted Operators, lift_operator.

Example

# Two qubits with ZZ coupling
sys1 = QuantumSystem([PAULIS[:X]], [(-1.0, 1.0)])
sys2 = QuantumSystem([PAULIS[:Y]], [(-1.0, 1.0)])
H_coupling = 0.1 * kron(PAULIS[:Z], PAULIS[:Z])
csys = CompositeQuantumSystem(H_coupling, [sys1, sys2], Float64[])

CompositeQuantumSystem(
    H_drift::AbstractMatrix,
    H_drives::AbstractVector{<:AbstractMatrix},
    subsystems::AbstractVector{<:QuantumSystem},
    drive_bounds::DriveBounds
)

Construct a CompositeQuantumSystem with coupling drift and drive terms.

Arguments

• H_drift::AbstractMatrix: Coupling drift Hamiltonian (in full tensor product space)
• H_drives::AbstractVector{<:AbstractMatrix}: Coupling drive Hamiltonians
• subsystems::AbstractVector{<:QuantumSystem}: Vector of subsystems to compose
• drive_bounds::DriveBounds: Drive bounds for the coupling drives (subsystem bounds are inherited). Can be:
  • Tuples (lower, upper) for asymmetric bounds
  • Scalars which are interpreted as symmetric bounds (-value, value)

The total drift includes both the coupling drift and all subsystem drifts (automatically lifted). The total drives include coupling drives followed by all subsystem drives (automatically lifted).

Example

sys1 = QuantumSystem(PAULIS[:Z], [PAULIS[:X]], [1.0])
sys2 = QuantumSystem([PAULIS[:Y]], [1.0])
g12 = 0.1 * kron(PAULIS[:X], PAULIS[:X])  # coupling drift
csys = CompositeQuantumSystem(g12, Matrix{ComplexF64}[], [sys1, sys2], Float64[])

CompositeQuantumSystem(
    subsystems::AbstractVector{<:QuantumSystem},
    drive_bounds::DriveBounds
)

Convenience constructor for a composite system with no coupling terms (neither drift nor drives).

Use this when you have independent subsystems that you want to represent in a single composite space, but without any direct coupling between them.

Arguments

• subsystems::AbstractVector{<:QuantumSystem}: Vector of subsystems to compose
• drive_bounds::DriveBounds: Drive bounds for the coupling drives (typically empty). Can be:
  • Tuples (lower, upper) for asymmetric bounds
  • Scalars which are interpreted as symmetric bounds (-value, value)

Example

sys1 = QuantumSystem([PAULIS[:X]], [1.0])
sys2 = QuantumSystem([PAULIS[:Y]], [1.0])
csys = CompositeQuantumSystem([sys1, sys2], Float64[])

CompositeQuantumSystem(
    H_drift::AbstractMatrix,
    subsystems::AbstractVector{<:QuantumSystem},
    drive_bounds::DriveBounds
)

Convenience constructor for a composite system with coupling drift but no coupling drives.

Arguments

• H_drift::AbstractMatrix: Coupling drift Hamiltonian
• subsystems::AbstractVector{<:QuantumSystem}: Vector of subsystems to compose
• drive_bounds::DriveBounds: Drive bounds for the coupling drives (typically empty). Can be:
  • Tuples (lower, upper) for asymmetric bounds
  • Scalars which are interpreted as symmetric bounds (-value, value)

Example

sys1 = QuantumSystem([PAULIS[:X]], [1.0])
sys2 = QuantumSystem([PAULIS[:Y]], [1.0])
H_coupling = 0.1 * kron(PAULIS[:Z], PAULIS[:Z])  # coupling drift
csys = CompositeQuantumSystem(H_coupling, [sys1, sys2], Float64[])

CompositeQuantumSystem(
    H_drives::AbstractVector{<:AbstractMatrix},
    subsystems::AbstractVector{<:QuantumSystem},
    drive_bounds::DriveBounds
)

Convenience constructor for a composite system with coupling drives but no coupling drift.

Arguments

• H_drives::AbstractVector{<:AbstractMatrix}: Coupling drive Hamiltonians
• subsystems::AbstractVector{<:QuantumSystem}: Vector of subsystems to compose
• drive_bounds::DriveBounds: Drive bounds for the coupling drives. Can be:
  • Tuples (lower, upper) for asymmetric bounds
  • Scalars which are interpreted as symmetric bounds (-value, value)

Example

sys1 = QuantumSystem([PAULIS[:X]], [1.0])
sys2 = QuantumSystem([PAULIS[:Y]], [1.0])
g12 = 0.1 * kron(PAULIS[:X], PAULIS[:X])  # coupling drive
csys = CompositeQuantumSystem([g12], [sys1, sys2], [1.0])  # symmetric bound

DriveBounds

Type alias for drive amplitude bounds input. Bounds can be specified as:

• A tuple (lower, upper) for asymmetric bounds
• A scalar value which is interpreted as symmetric bounds (-value, value)

Examples

drive_bounds = [(-1.0, 1.0), 0.5, (-0.3, 0.7)]
# Interpreted as: [(-1.0, 1.0), (-0.5, 0.5), (-0.3, 0.7)]

LinearDrive <: AbstractDrive

Standard linear drive: coefficient is u[index].

This is the default representation when constructing QuantumSystem(H_drift, H_drives, bounds).

Fields

• H::SparseMatrixCSC{ComplexF64,Int}: The Hermitian drive operator
• index::Int: Index into the control vector u

Example

LinearDrive(sparse(PAULIS.X), 1)  # u[1] * σx

NonlinearDrive{F,DF} <: AbstractDrive

Drive term with a nonlinear scalar coefficient: f(u) * H.

The user provides either:

1. Both the coefficient function and its Jacobian (3-arg form), or
2. Just the coefficient function (2-arg form) — the Jacobian is computed automatically via ForwardDiff.

Fields

• H::SparseMatrixCSC{ComplexF64,Int}: The Hermitian drive operator
• coeff::F: (u::AbstractVector) -> scalar — coefficient function
• coeff_jac::DF: (u::AbstractVector, j::Int) -> scalar — ∂coeff/∂u_j
• active_controls::Vector{Int}: Control indices where ∂coeff/∂u_j can be nonzero. Empty means "all controls" (no structural sparsity info).

Example: Auto-Jacobian (recommended)

NonlinearDrive(σz / 2, u -> u[3]^2 + u[4]^2)

Example: Manual Jacobian with active controls

NonlinearDrive(
    σz / 2,
    u -> u[3]^2 + u[4]^2,
    (u, j) -> j == 3 ? 2u[3] : j == 4 ? 2u[4] : 0.0;
    active_controls = [3, 4]
)

NonlinearDrive(H::AbstractMatrix, coeff, coeff_jac; active_controls=Int[])

Construct a NonlinearDrive with an explicit Jacobian function. active_controls lists which control indices have nonzero ∂coeff/∂u_j (empty = all).

NonlinearDrive(H::AbstractMatrix, coeff; active_controls=Int[])

Construct a NonlinearDrive with an auto-generated Jacobian via ForwardDiff.

This is the recommended constructor — the Jacobian is computed automatically from the coefficient function, eliminating the risk of hand-written Jacobian errors.

Example

# Displaced-frame Stark shift: coefficient = u₃² + u₄²
drive = NonlinearDrive(σz / 2, u -> u[3]^2 + u[4]^2; active_controls = [3, 4])

OpenQuantumSystem <: AbstractQuantumSystem

A struct for storing open quantum dynamics.

Fields

• H::Function: The Hamiltonian function: (u, t) -> H(u, t)
• 𝒢::Function: The Lindbladian generator function: u -> 𝒢(u)
• H_drift::SparseMatrixCSC{ComplexF64, Int}: The drift Hamiltonian
• H_drives::Vector{AbstractDrive}: The canonical drive terms, each pairing an operator with a coefficient function. Matrix-based constructors auto-populate this with LinearDrive objects; function-based systems leave it empty.
• drive_bounds::Vector{Tuple{Float64, Float64}}: Drive amplitude bounds
• n_drives::Int: The number of control drives
• levels::Int: The number of levels in the system
• dissipation_operators::Vector{SparseMatrixCSC{ComplexF64, Int}}: The dissipation operators
• time_dependent::Bool: Whether the Hamiltonian has explicit time dependence
• global_params::NamedTuple: Global parameters stored with the system (e.g., physical constants)

See also QuantumSystem.

OpenQuantumSystem(
    H_drift::AbstractMatrix,
    drives::Vector{<:AbstractDrive},
    drive_bounds::DriveBounds;
    dissipation_operators=Matrix{ComplexF64}[],
    time_dependent::Bool=false,
    global_params::NamedTuple=NamedTuple()
)

Construct an OpenQuantumSystem from a drift Hamiltonian and typed drive terms.

This constructor supports both linear and nonlinear drives. The control dimension is determined by length(drive_bounds), which may differ from length(drives) when nonlinear drives combine multiple controls.

Example

L_ops = [sqrt(γ) * annihilate(levels)]
drives = [
    LinearDrive(sparse(σx), 1),
    NonlinearDrive(σz, u -> u[1]^2 + u[2]^2),
]
sys = OpenQuantumSystem(H_drift, drives, [1.0, 1.0]; dissipation_operators=L_ops)

OpenQuantumSystem(
    H_drift::AbstractMatrix{<:Number},
    H_drives::AbstractVector{<:AbstractMatrix{<:Number}},
    drive_bounds::DriveBounds;
    dissipation_operators=Matrix{ComplexF64}[],
    time_dependent::Bool=false,
    global_params::NamedTuple=NamedTuple()
)
OpenQuantumSystem(
    H_drift::AbstractMatrix{<:Number};
    dissipation_operators=Matrix{ComplexF64}[],
    time_dependent::Bool=false,
    global_params::NamedTuple=NamedTuple()
)
OpenQuantumSystem(
    H_drives::Vector{<:AbstractMatrix{<:Number}},
    drive_bounds::DriveBounds;
    dissipation_operators=Matrix{ComplexF64}[],
    time_dependent::Bool=false,
    global_params::NamedTuple=NamedTuple()
)
OpenQuantumSystem(
    H::Function,
    drive_bounds::DriveBounds;
    dissipation_operators=Matrix{ComplexF64}[],
    time_dependent::Bool=false,
    global_params::NamedTuple=NamedTuple()
)
OpenQuantumSystem(
    system::QuantumSystem;
    dissipation_operators=Matrix{ComplexF64}[]
)

Constructs an OpenQuantumSystem object from the drift and drive Hamiltonian terms and dissipation operators.

Drive Bounds

The drive_bounds parameter can be:

• Tuples (lower, upper) for asymmetric bounds
• Scalars which are interpreted as symmetric bounds (-value, value)

OpenQuantumSystem(system::QuantumSystem; dissipation_operators=[])

Construct an OpenQuantumSystem from a QuantumSystem by adding dissipation operators.

When the QuantumSystem has typed drives (including nonlinear), they are preserved in the resulting OpenQuantumSystem.

OpenQuantumSystem(H_drift; dissipation_operators=[], time_dependent=false)

Construct an OpenQuantumSystem with only drift (no drives).

OpenQuantumSystem(H::Function, drive_bounds; dissipation_operators=[], time_dependent=false)

Construct an OpenQuantumSystem from a Hamiltonian function.

OpenQuantumSystem(H_drives, drive_bounds; dissipation_operators=[], time_dependent=false)

Construct an OpenQuantumSystem with no drift.

QuantumSystem <: AbstractQuantumSystem

A struct for storing quantum dynamics.

Fields

• H::Function: The Hamiltonian function: (u, t) -> H(u, t), where u is the control vector and t is time
• G::Function: The isomorphic generator function: (u, t) -> G(u, t), including the Hamiltonian mapped to superoperator space
• H_drift::SparseMatrixCSC{ComplexF64, Int}: The drift Hamiltonian (time-independent component)
• H_drives::Vector{AbstractDrive}: The canonical drive terms, each pairing an operator with a coefficient function. Matrix-based constructors auto-populate this with LinearDrive objects; function-based systems leave it empty.
• drive_bounds::Vector{Tuple{Float64, Float64}}: Drive amplitude bounds for each control (lower, upper)
• n_drives::Int: The number of control channels (length of the control vector u)
• levels::Int: The number of levels (dimension) in the system
• time_dependent::Bool: Whether the Hamiltonian has explicit time dependence beyond control modulation
• global_params::NamedTuple: Global parameters that the Hamiltonian may depend on (e.g., (δ=0.5, Ω=1.0))

See also OpenQuantumSystem, VariationalQuantumSystem.

QuantumSystem(
    H_drift::AbstractMatrix,
    drives::Vector{<:AbstractDrive},
    drive_bounds::Vector;
    time_dependent::Bool=false,
    global_params::NamedTuple=NamedTuple()
)

Construct a QuantumSystem from a drift Hamiltonian and typed drive terms.

This constructor supports both linear and nonlinear drives. The control dimension is determined by length(drive_bounds), which may differ from length(drives) when nonlinear drives combine multiple controls.

The resulting Hamiltonian is: H(u, t) = Hdrift + Σd drive_coeff(d, u) * d.H

Arguments

• H_drift::AbstractMatrix: The drift (time-independent) Hamiltonian
• drives::Vector{<:AbstractDrive}: Vector of drive terms (LinearDrive or NonlinearDrive)
• drive_bounds::DriveBounds: Bounds for each physical control. Length = control dimension.

Example: Displaced frame with nonlinear |α|² term

sys = QuantumSystem(
    H_drift,
    [
        LinearDrive(sparse(σx), 1),                    # u[1] * σx (qubit I)
        LinearDrive(sparse(σy), 2),                    # u[2] * σy (qubit Q)
        LinearDrive(sparse(χ * Xa * σz), 3),           # u[3] * χ·Xa·σz (displacement I)
        LinearDrive(sparse(χ * Pa * σz), 4),           # u[4] * χ·Pa·σz (displacement Q)
        NonlinearDrive(                                 # (u[3]²+u[4]²) * χ·σz/2
            sparse(χ * σz / 2),
            u -> u[3]^2 + u[4]^2,
            (u, j) -> j == 3 ? 2u[3] : j == 4 ? 2u[4] : 0.0
        ),
    ],
    [Ω_bound, Ω_bound, α_bound, α_bound]  # 4 physical controls
)

QuantumSystem(
    H_drift::AbstractMatrix{<:Number},
    H_drives::Vector{<:AbstractMatrix{<:Number}},
    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}};
    time_dependent::Bool=false
)

Construct a QuantumSystem from drift and drive Hamiltonian terms.

Arguments

• H_drift::AbstractMatrix: The drift (time-independent) Hamiltonian
• H_drives::Vector{<:AbstractMatrix}: Vector of drive Hamiltonians, one for each control
• drive_bounds::DriveBounds: Drive amplitude bounds for each control. Can be:
  • Tuples (lower, upper) for asymmetric bounds
  • Scalars which are interpreted as symmetric bounds (-value, value)

Keyword Arguments

• time_dependent::Bool=false: Set to true if using time-dependent modulation (typically handled at a higher level)
• global_params::NamedTuple=NamedTuple(): Global parameters stored with the system. Note: for matrix-based systems, matrices are fixed at construction, so globalparams are mainly for storage/bookkeeping and later updates via `updateglobal_params!`

The resulting Hamiltonian is: H(u, t) = Hdrift + Σᵢ uᵢ * Hdrives[i]

Example

sys = QuantumSystem(
    PAULIS[:Z],                    # drift
    [PAULIS[:X], PAULIS[:Y]],      # drives
    [1.0, 1.0]                     # symmetric bounds: [(-1.0, 1.0), (-1.0, 1.0)]
)

QuantumSystem(H::Function, drive_bounds::Vector; time_dependent::Bool=false)

Construct a QuantumSystem from a Hamiltonian function.

Arguments

• H::Function: Hamiltonian function with signature (u, t) -> H(u, t) where:
  • u is a vector containing [controls..., globals...] (if system has global parameters)
  • For matrix-based systems, only the first n_drives elements are used for controls
  • For function-based systems, handle globals via closure or by accessing u beyond control indices
  • t is time
• drive_bounds::DriveBounds: Drive amplitude bounds for each control. Can be:
  • Tuples (lower, upper) for asymmetric bounds
  • Scalars which are interpreted as symmetric bounds (-value, value)

Keyword Arguments

• time_dependent::Bool=false: Set to true if the Hamiltonian has explicit time dependence (e.g., cos(ωt) modulation)
• global_params::NamedTuple=NamedTuple(): Global parameters stored with the system for bookkeeping

Example

# Define a time-dependent Hamiltonian
H = (u, t) -> PAULIS[:Z] + u[1] * cos(ω * t) * PAULIS[:X]
sys = QuantumSystem(H, [(-1.0, 1.0)]; time_dependent=true)

QuantumSystem(drives::Vector{<:AbstractDrive}, drive_bounds::Vector; time_dependent::Bool=false)

Convenience constructor for a typed-drive system with no drift Hamiltonian (H_drift = 0).

Example

sys = QuantumSystem([LinearDrive(sparse(PAULIS[:X]), 1)], [1.0])

QuantumSystem(H_drift::AbstractMatrix; time_dependent::Bool=false)

Convenience constructor for a system with only a drift Hamiltonian (no drives).

Example

sys = QuantumSystem(PAULIS[:Z])

QuantumSystem(H_drives::Vector{<:AbstractMatrix}, drive_bounds::Vector; time_dependent::Bool=false)

Convenience constructor for a system with no drift Hamiltonian (H_drift = 0).

Arguments

• H_drives::Vector{<:AbstractMatrix}: Vector of drive Hamiltonians
• drive_bounds::DriveBounds: Drive amplitude bounds for each control. Can be:
  • Tuples (lower, upper) for asymmetric bounds
  • Scalars which are interpreted as symmetric bounds (-value, value)

Example

# Using scalars for symmetric bounds
sys = QuantumSystem([PAULIS[:X], PAULIS[:Y]], [1.0, 1.0])
# Equivalent to: drive_bounds = [(-1.0, 1.0), (-1.0, 1.0)]

VariationalQuantumSystem <: AbstractQuantumSystem

A struct for storing variational quantum dynamics, used for sensitivity and robustness analysis.

Variational systems allow exploring how the dynamics change under perturbations to the Hamiltonian. The variational operators represent directions of uncertainty or perturbation in the system.

Fields

• H::Function: The Hamiltonian function: (u, t) -> H(u, t)
• G::Function: The isomorphic generator function: (u, t) -> G(u, t)
• G_vars::AbstractVector{<:Function}: Variational generator functions, one for each perturbation direction
• drive_bounds::Vector{Tuple{Float64, Float64}}: Drive amplitude bounds
• n_drives::Int: The number of control drives in the system
• levels::Int: The number of levels (dimension) in the system

See also QuantumSystem, OpenQuantumSystem.

VariationalQuantumSystem(
    H_drift::AbstractMatrix,
    H_drives::AbstractVector{<:AbstractMatrix},
    H_vars::AbstractVector{<:AbstractMatrix},
    drive_bounds::DriveBounds
)

Construct a VariationalQuantumSystem from drift, drive, and variational Hamiltonian terms.

Arguments

• H_drift::AbstractMatrix: The drift (time-independent) Hamiltonian
• H_drives::AbstractVector{<:AbstractMatrix}: Vector of drive Hamiltonians for control
• H_vars::AbstractVector{<:AbstractMatrix}: Vector of variational Hamiltonians representing perturbation directions
• drive_bounds::DriveBounds: Drive amplitude bounds for each control. Can be:
  • Tuples (lower, upper) for asymmetric bounds
  • Scalars which are interpreted as symmetric bounds (-value, value)

The variational operators allow sensitivity analysis by exploring how dynamics change under perturbations: Hperturbed = H + Σᵢ εᵢ * Hvars[i]

Example

varsys = VariationalQuantumSystem(
    PAULIS[:Z],                    # drift
    [PAULIS[:X], PAULIS[:Y]],      # drives
    [PAULIS[:X]],                  # variational perturbations
    [1.0, 1.0]                     # symmetric bounds: [(-1.0, 1.0), (-1.0, 1.0)]
)

Isomorphisms.G(d::AbstractDrive)

Delegate G to the underlying Hamiltonian matrix d.H, so that broadcasting G.(sys.H_drives) works when H_drives contains AbstractDrive objects.

_float_params(nt::NamedTuple)

Convert all values in a NamedTuple to their floating-point equivalents. Ensures type-stable ODE solutions when global parameters are updated during optimization.

_forwarddiff_jacobian(f) -> (u, j) -> ∂f/∂u_j

Create a Jacobian function from a scalar-valued function f(u) using ForwardDiff.

active_controls(d::AbstractDrive) -> Vector{Int}

Return the control indices where drive_coeff_jac(d, u, j) can be nonzero. An empty vector means "all controls" (no structural sparsity information available).

For LinearDrive, returns [d.index]. For NonlinearDrive, returns d.active_controls (user-specified or empty).

compact_lindbladian_generators(sys::OpenQuantumSystem)

Compute the compact Lindbladian generators for use with the compact density isomorphism. Returns (𝒢c_drift, 𝒢c_drives) where:

• 𝒢c_drift = P * (𝒢_drift + 𝒟) * L — compact drift + dissipation generator (n² × n²)
• 𝒢c_drives[i] — compact generator for each drive term (each n² × n²)

For linear-only systems, 𝒢c_drives[i] corresponds to control u[i]. For systems with nonlinear drives, 𝒢c_drives[i] corresponds to drive term i, and must be weighted by drive_coeff(sys.H_drives[i], u) instead of u[i].

drive_coeff(d::AbstractDrive, u::AbstractVector) -> Number

Compute the scalar coefficient of this drive at controls u.

drive_coeff_jac(d::AbstractDrive, u::AbstractVector, j::Int) -> Number

Compute ∂coeff/∂u_j at controls u.

For LinearDrive: returns 1.0 if j == d.index, 0.0 otherwise (Kronecker delta). For NonlinearDrive: evaluates the Jacobian function (user-provided or auto-generated).

get_drift(sys::AbstractQuantumSystem)

Returns the drift Hamiltonian of the system.

get_drive_terms(sys::AbstractQuantumSystem) -> Vector{AbstractDrive}

Return the typed drive terms from the system. Each AbstractDrive pairs a Hamiltonian operator with a scalar coefficient function and Jacobian.

Returns sys.H_drives directly when available, or an empty vector for function-based systems that don't use typed drives.

See also get_drives for just the Hamiltonian matrices.

get_drives(sys::AbstractQuantumSystem)

Returns the drive Hamiltonian matrices of the system.

For systems with typed drives (sys.H_drives), returns the H matrix from each drive term. For function-based systems without drives, extracts operators via basis vector evaluation.

has_nonlinear_drives(drives::Vector{<:AbstractDrive}) -> Bool

Check if any drive terms are nonlinear.

is_hermitian(H::AbstractMatrix; tol=1e-10)

Check if a matrix is Hermitian within a tolerance.

normalize_drive_bounds(bounds::DriveBounds)

Convert drive bounds to a consistent tuple format. Scalar values are converted to symmetric bounds around zero: b becomes (-b, b).

Arguments

• bounds::DriveBounds: Input bounds, can be tuples or scalars

Returns

• Vector{Tuple{Float64, Float64}}: Normalized bounds as tuples

Examples

# All scalars (symmetric bounds)
normalize_drive_bounds([1.0, 1.5, 0.5])
# Returns: [(-1.0, 1.0), (-1.5, 1.5), (-0.5, 0.5)]

# All tuples (asymmetric bounds)
normalize_drive_bounds([(-2.0, 3.0), (-1.0, 1.0)])
# Returns: [(-2.0, 3.0), (-1.0, 1.0)]

# Mixed types (requires explicit type annotation)
normalize_drive_bounds(Union{Float64, Tuple{Float64,Float64}}[1.0, (-2.0, 3.0)])
# Returns: [(-1.0, 1.0), (-2.0, 3.0)]

validate_drive_jacobian(d::NonlinearDrive, n_controls::Int; atol=1e-6, n_samples=3)

Spot-check the Jacobian of a NonlinearDrive against ForwardDiff at random control vectors. Throws an AssertionError if the user-provided Jacobian disagrees with the AD Jacobian.

This is called automatically during QuantumSystem construction for all NonlinearDrive terms, catching sign errors or off-by-one bugs early.

CatSystem(;
    g2=0.36, χ_aa=-7e-3, χ_bb=-32, χ_ab=0.79,
    κa=53e-3, κb=13,
    cat_levels=13, buffer_levels=3,
    prefactor=1, drive_bounds=[1.0, 1.0],
)::OpenQuantumSystem

Construct an OpenQuantumSystem for a two-mode cat qubit (cat ⊗ buffer).

Hamiltonian

The drift Hamiltonian includes Kerr, cross-Kerr, and two-photon exchange terms:

\[H = -\frac{\chi_{aa}}{2} a^{\dagger 2} a^2 -\frac{\chi_{bb}}{2} b^{\dagger 2} b^2 -\chi_{ab}\, a^\dagger a\, b^\dagger b + g_2\, a^{\dagger 2} b + g_2^*\, a^2 b^\dagger\]

The two drives are the buffer displacement $b + b^\dagger$ and a Kerr correction $a^\dagger a$.

Keyword Arguments

• g2: Two-photon exchange coupling (MHz · 2π)
• χ_aa: Cat self-Kerr (MHz · 2π)
• χ_bb: Buffer self-Kerr (MHz · 2π)
• χ_ab: Cross-Kerr between cat and buffer (MHz · 2π)
• κa: Cat decay rate (MHz · 2π)
• κb: Buffer decay rate (MHz · 2π)
• cat_levels: Truncation of cat mode Fock space
• buffer_levels: Truncation of buffer mode Fock space
• prefactor: Global scaling applied to all couplings and rates
• drive_bounds: Bounds on the two drive amplitudes

All parameters are scaled by $2π$ (Hamiltonian) or $\\sqrt{2π}$ (dissipators) internally.

IonChainSystem(;
    N_ions::Int=2,
    ion_levels::Int=2,
    N_modes::Int=1,
    mode_levels::Int=10,
    ωq::Union{Float64, Vector{Float64}}=1.0,
    ωm::Union{Float64, Vector{Float64}}=0.1,
    η::Union{Float64, Matrix{Float64}}=0.1,
    lab_frame::Bool=false,
    frame_ω::Float64=lab_frame ? 0.0 : (ωq isa Vector ? ωq[1] : ωq),
    multiply_by_2π::Bool=true,
    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}=fill(1.0, 2*N_ions),
) -> QuantumSystem

Returns a QuantumSystem object for a chain of trapped ions coupled via motional modes.

The system consists of N_ions ions, each with ion_levels internal states, coupled to N_modes shared motional modes with mode_levels Fock states each.

Hamiltonian

In the lab frame:

\[H = \sum_{i=1}^{N_{\text{ions}}} \omega_{q,i} \sigma_i^+ \sigma_i^- + \sum_{m=1}^{N_{\text{modes}}} \omega_{m} a_m^\dagger a_m + \sum_{i,m} \eta_{i,m} (\sigma_i^+ + \sigma_i^-)(a_m + a_m^\dagger) + \sum_i \Omega_{x,i}(t) \sigma_i^x + \sum_i \Omega_{y,i}(t) \sigma_i^y\]

In the rotating frame at frequency frame_ω:

\[H = \sum_{i=1}^{N_{\text{ions}}} (\omega_{q,i} - \omega_{\text{frame}}) \sigma_i^+ \sigma_i^- + \sum_{m=1}^{N_{\text{modes}}} \omega_{m} a_m^\dagger a_m + \sum_{i,m} \eta_{i,m} (\sigma_i^+ + \sigma_i^-)(a_m + a_m^\dagger) + \sum_i \Omega_{x,i}(t) \sigma_i^x + \sum_i \Omega_{y,i}(t) \sigma_i^y\]

where:

• $\sigma_i^+, \sigma_i^-$ are raising/lowering operators for ion i
• $\sigma_i^x, \sigma_i^y$ are Pauli operators for ion i
• $a_m, a_m^\dagger$ are annihilation/creation operators for mode m
• $\omega_{q,i}$ is the transition frequency of ion i
• $\omega_m$ is the frequency of motional mode m
• $\eta_{i,m}$ is the Lamb-Dicke parameter coupling ion i to mode m

Keyword Arguments

• N_ions: Number of ions in the chain
• ion_levels: Number of internal levels per ion (typically 2 for qubit)
• N_modes: Number of motional modes to include
• mode_levels: Number of Fock states per motional mode
• ωq: Ion transition frequency (or frequencies). Scalar or vector of length N_ions. In GHz.
• ωm: Motional mode frequency (or frequencies). Scalar or vector of length N_modes. In GHz.
• η: Lamb-Dicke parameter(s). Scalar (uniform coupling), or N_ions × N_modes matrix.
• lab_frame: If true, use lab frame Hamiltonian. If false, use rotating frame.
• frame_ω: Rotating frame frequency in GHz. Defaults to first ion frequency.
• multiply_by_2π: Whether to multiply Hamiltonian by 2π (default true, since frequencies are in GHz).
• drive_bounds: Control bounds. Vector of length 2*N_ions for [Ωx₁, Ωy₁, Ωx₂, Ωy₂, ...].

Example

# Two ions, one motional mode, Mølmer-Sørensen setup
sys = IonChainSystem(
    N_ions=2,
    N_modes=1,
    ωq=1.0,      # 1 GHz qubit frequency
    ωm=0.1,      # 100 MHz mode frequency  
    η=0.1,       # Lamb-Dicke parameter
    mode_levels=5,
)

References

• Sørensen, A. & Mølmer, K. "Quantum computation with ions in thermal motion." Phys. Rev. Lett. 82, 1971 (1999).
• Sørensen, A. & Mølmer, K. "Entanglement and quantum computation with ions in thermal motion." Phys. Rev. A 62, 022311 (2000).

MolmerSorensenCoupling(
    N_ions::Int,
    N_modes::Int,
    ion_levels::Int,
    mode_levels::Int,
    η::Union{Float64, Matrix{Float64}},
    ωm::Union{Float64, Vector{Float64}},
) -> Matrix{ComplexF64}

Returns the Mølmer-Sørensen coupling term for an ion chain, which mediates effective ion-ion interactions via the motional modes.

In the Lamb-Dicke regime with appropriate drive detunings, the effective Hamiltonian is:

\[H_{\text{MS}} = \sum_{i<j} J_{ij} \sigma_i^x \sigma_j^x\]

where the coupling strength is:

\[J_{ij} = \sum_m \frac{\eta_{i,m} \eta_{j,m} \Omega_i \Omega_j}{4 \Delta_m}\]

with $\Delta_m$ being the detuning from mode $m$.

Arguments

• N_ions: Number of ions
• N_modes: Number of motional modes
• ion_levels: Internal levels per ion
• mode_levels: Fock states per mode
• η: Lamb-Dicke parameters (scalar or Nions × Nmodes matrix)
• ωm: Mode frequencies (scalar or vector)

Returns

Matrix representing the σˣᵢ σˣⱼ interaction for use in building MS gates.

MultiTransmonSystem(
    ωs::AbstractVector{Float64},
    δs::AbstractVector{Float64},
    gs::AbstractMatrix{Float64};
    drive_bounds::Union{Float64, Vector{<:Union{Tuple{Float64, Float64}, Float64}}}=1.0,
    levels_per_transmon::Int = 3,
    subsystem_levels::AbstractVector{Int} = fill(levels_per_transmon, length(ωs)),
    lab_frame=false,
    subsystems::AbstractVector{Int} = 1:length(ωs),
    subsystem_drive_indices::AbstractVector{Int} = 1:length(ωs),
    kwargs...
) -> CompositeQuantumSystem

Returns a CompositeQuantumSystem object for a multi-transmon system.

RadialMSGateSystem(;
    N_ions::Int=2,
    mode_levels::Int=5,
    ωm_radial::Vector{Float64}=[5.0, 5.0, 5.1, 5.1],  # 4 radial modes for 2 ions
    δ::Union{Float64, Vector{Float64}}=0.2,           # Detuning(s) from mode(s)
    η::Union{Float64, Matrix{Float64}}=0.1,           # Lamb-Dicke parameters
    multiply_by_2π::Bool=true,
    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}=fill(1.0, N_ions),
) -> QuantumSystem

Returns a time-dependent QuantumSystem for the radial-mode Mølmer-Sørensen gate as described in the paper:

"Realization and Calibration of Continuously Parameterized Two-Qubit Gates on a Trapped-Ion Quantum Processor" (IEEE TQE 2024)

This implements the MS gate using only radial motional modes (not axial modes), which provides 2N modes for N ions (N modes along each of two transverse axes).

Hamiltonian (Equation 2 from paper)

In the interaction picture:

\[H(t) = -\frac{i\hbar}{2} \sum_{i,k} \sigma_{x,i} \eta_{k,i} \Omega_i a_k e^{-i\delta_k t} + \text{h.c.}\]

Expanding the Hermitian conjugate:

\[H(t) = -\frac{i}{2} \sum_{i,k} \eta_{k,i} \Omega_i \sigma_{x,i} \left(a_k e^{-i\delta_k t} - a_k^\dagger e^{i\delta_k t}\right)\]

where:

• \[k\]

  indexes the 2N radial modes (N along x, N along y)

• \[\sigma_{x,i}\]

  is Pauli-X on ion $i$

• \[\eta_{k,i}\]

  is the Lamb-Dicke parameter for ion $i$, mode $k$

• \[\Omega_i(t)\]

  is the control amplitude (Rabi frequency) for ion $i$

• \[\delta_k\]

  is the detuning from motional sideband of mode $k$

• \[a_k, a_k^\dagger\]

  are phonon operators for radial mode $k$

Radial Mode Structure

For N ions in a linear trap with radial confinement:

• Axial modes (along trap axis): Not used for this gate
• Radial modes: 2N modes total
  • N modes along transverse x-direction
  • N modes along transverse y-direction

For N=2 ions: 4 radial modes participate in the gate dynamics.

Typical Parameters (Q-SCOUT platform at Sandia, ¹⁷¹Yb⁺)

• Radial frequencies: $\omega_r / 2\pi \sim 5$ MHz (higher than axial ~2 MHz)
• Lamb-Dicke: $\eta \sim 0.05 - 0.15$
• Detuning: $\delta / 2\pi \sim 100 - 500$ kHz
• Gate time: $50 - 200$ μs
• Phonon states: $n_{\max} = 3-5$ typically sufficient

Keyword Arguments

• N_ions: Number of ions (default: 2)
• mode_levels: Fock states per radial mode (default: 5)
• ωm_radial: Radial mode frequencies in GHz. Vector of length 2N. Example for 2 ions: [5.0, 5.0, 5.1, 5.1] (nearly degenerate pairs)
• δ: Detuning(s) from sideband in GHz. Scalar (uniform) or vector per mode.
• η: Lamb-Dicke parameter(s). Scalar (uniform) or N_ions × 2N matrix.
• multiply_by_2π: Multiply by 2π (default true, since frequencies in GHz)
• drive_bounds: Control amplitude bounds for each ion (length N_ions)

Example: Two-Ion Radial MS Gate

sys = RadialMSGateSystem(
    N_ions=2,
    mode_levels=5,
    ωm_radial=[5.0, 5.0, 5.1, 5.1],  # Two nearly-degenerate pairs
    δ=0.2,                            # 200 kHz detuning
    η=0.1,                            # Lamb-Dicke parameter
    drive_bounds=[1.0, 1.0]           # Amplitude bounds (GHz)
)

# Create trajectory for XX gate
U_goal = exp(-im * π/4 * kron([0 1; 1 0], [0 1; 1 0]))  # XX(π/2)
qtraj = UnitaryTrajectory(sys, U_goal, 100e-6)  # 100 μs gate

Optimization Considerations

1. Motional closure: All 2N modes must satisfy $|\alpha_k(\tau)| \approx 0$
2. Target mode: Choose one mode (e.g., k=1) as primary entangling mode
3. Spectator modes: Other modes should remain minimally excited
4. Control strategy: Individual ion addressing via $\Omega_i(t)$

References

• Sørensen & Mølmer, "Quantum computation with ions in thermal motion," PRL 82, 1971 (1999)
• Mizrahi et al., "Realization and Calibration of Continuously Parameterized Two-Qubit Gates...," IEEE TQE (2024)

RadialMSGateSystemWithPhase(;
    N_ions::Int=2,
    mode_levels::Int=5,
    ωm_radial::Vector{Float64}=[5.0, 5.0, 5.1, 5.1],
    δ::Union{Float64, Vector{Float64}}=0.2,
    η::Union{Float64, Matrix{Float64}}=0.1,
    multiply_by_2π::Bool=true,
    amplitude_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}=fill(1.0, N_ions),
    phase_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}=fill((-π, π), N_ions),
) -> QuantumSystem

Returns a time-dependent QuantumSystem for the radial-mode MS gate with phase controls to enable AC Stark shift compensation.

Hamiltonian (with phase modulation)

Instead of $\sigma_x = \sigma^+ + \sigma^-$, we use phase-modulated drives:

\[H(t) = \frac{1}{2} \sum_{i,k} \eta_{k,i} \Omega_i(t) \left(\sigma^+_i e^{i\phi_i(t)} + \sigma^-_i e^{-i\phi_i(t)}\right) \left(a_k e^{-i\delta_k t} + a_k^\dagger e^{i\delta_k t}\right)\]

where $\Omega_i(t)$ and $\phi_i(t)$ are independent controls.

Why Phase Controls?

Off-resonant coupling to spectator modes creates AC Stark shifts:

\[\Delta E_{\text{Stark}} \sim \frac{\eta^2 \Omega^2(t)}{\delta_{\text{spectator}}}\]

This causes time-varying phase accumulation that $\sigma_x$ control alone cannot compensate. The solution: actively modulate $\phi_i(t)$ to cancel the Stark-induced phase, typically using:

\[\phi(t) \sim \int_0^t \frac{\eta^2 \Omega^2(t')}{\delta} dt' \sim \text{erf}(\sqrt{2}t) \text{ for Gaussian pulses}\]

This enables loop closure in phase space and high-fidelity gates ($F > 0.99$).

Control Structure

Controls: $[\\Omega_1, \phi_1, \\Omega_2, \phi_2, \ldots]$ for $N_{\text{ions}}$ ions.

• Even indices (1, 3, 5, ...): Amplitudes $\Omega_i(t)$ (Rabi frequency)
• Odd indices (2, 4, 6, ...): Phases $\phi_i(t)$ (laser phase)

Example

sys = RadialMSGateSystemWithPhase(
    N_ions=2,
    mode_levels=3,
    ωm_radial=[5.0, 5.0, 5.1, 5.1],
    δ=0.2,
    η=0.1,
    amplitude_bounds=[1.0, 1.0],
    phase_bounds=[(-Float64(π), Float64(π)), (-Float64(π), Float64(π))]
)

# sys.n_drives == 4 (Ω₁, φ₁, Ω₂, φ₂)

See Also

• RadialMSGateSystem: Amplitude-only version (simpler but limited fidelity)

RydbergChainSystem(;
    N::Int=3, # number of atoms
    C::Float64=862690*2π,
    distance::Float64=10.0, # μm
    cutoff_order::Int=2, # 1 is nearest neighbor, 2 is next-nearest neighbor, etc.
    local_detune::Bool=false, # If true, include one local detuning pattern.
    all2all::Bool=true, # If true, include all-to-all interactions.
    ignore_Y_drive::Bool=false, # If true, ignore the Y drive. (In the experiments, X&Y drives are implemented by Rabi amplitude and its phase.)
    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}=[1.0, 1.0, 1.0], # Bounds for [Ωx, Ωy, Δ] or [Ωx, Δ] if ignore_Y_drive
) -> QuantumSystem

Returns a QuantumSystem object for the Rydberg atom chain in the spin basis |g⟩ = |0⟩ = [1, 0], |r⟩ = |1⟩ = [0, 1].

\[H = \sum_i 0.5*\Omega_i(t)\cos(\phi_i(t)) \sigma_i^x - 0.5*\Omega_i(t)\sin(\phi_i(t)) \sigma_i^y - \sum_i \Delta_i(t)n_i + \sum_{i<j} \frac{C}{|i-j|^6} n_i n_j\]

Keyword Arguments

• N: Number of atoms.
• C: The Rydberg interaction strength in MHz*μm^6.
• distance: The distance between atoms in μm.
• cutoff_order: Interaction range cutoff, 1 is nearest neighbor, 2 is next nearest neighbor.
• local_detune: If true, include one local detuning pattern.
• all2all: If true, include all-to-all interactions.
• ignore_Y_drive: If true, ignore the Y drive. (In the experiments, X&Y drives are implemented by Rabi amplitude and its phase.)
• drive_bounds: Bounds for drive amplitudes.

TransmonCavitySystem(;
    qubit_levels::Int=4,
    cavity_levels::Int=12,
    χ::Float64=2π * 32.8e-6,    # Dispersive shift (GHz)
    χ′::Float64=2π * 1.5e-9,    # Higher-order dispersive shift (GHz)
    K_c::Float64=2π * 1e-9 / 2, # Cavity self-Kerr (GHz)
    K_q::Float64=2π * 193e-3 / 2, # Qubit self-Kerr (GHz)
    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}=fill(1.0, 4),
    multiply_by_2π::Bool=false, # Already in GHz with 2π factors
) -> QuantumSystem

Returns a QuantumSystem for a transmon qubit dispersively coupled to a cavity mode.

This system models circuit QED architectures where a transmon (artificial atom) is coupled to a microwave resonator in the dispersive regime, enabling quantum state manipulation and readout.

Hamiltonian

\[H = \tilde{\Delta} \hat{b}^\dagger \hat{b} - \chi \hat{a}^\dagger \hat{a} \hat{b}^\dagger \hat{b} - \chi' \hat{b}^{\dagger 2} \hat{b}^2 \hat{a}^\dagger \hat{a} - K_q \hat{a}^{\dagger 2} \hat{a}^2 - K_c \hat{b}^{\dagger 2} \hat{b}^2\]

where:

• $\hat{a}$, $\hat{a}^\dagger$ are the transmon annihilation/creation operators
• $\hat{b}$, $\hat{b}^\dagger$ are the cavity annihilation/creation operators
• $\tilde{\Delta} = \chi/2$ is the shifted cavity frequency
• $\chi$ is the dispersive shift (qubit-cavity coupling)
• $\chi'$ is a higher-order dispersive correction
• $K_q$, $K_c$ are self-Kerr nonlinearities

The drives are:

1. $\hat{a}^\dagger + \hat{a}$ - Real transmon drive
2. $i(\hat{a}^\dagger - \hat{a})$ - Imaginary transmon drive
3. $\hat{b}^\dagger + \hat{b}$ - Real cavity drive
4. $i(\hat{b}^\dagger - \hat{b})$ - Imaginary cavity drive

Keyword Arguments

• qubit_levels: Number of transmon Fock states (typically 3-5)
• cavity_levels: Number of cavity Fock states (typically 10-20)
• χ: Dispersive shift in GHz (with 2π). Typical: ~2π × 30-50 kHz
• χ′: Higher-order dispersive shift in GHz. Typically small (~2π × 1-2 Hz)
• K_c: Cavity self-Kerr in GHz. Typically ~2π × 1 Hz
• K_q: Qubit self-Kerr (anharmonicity/2) in GHz. Typical: ~2π × 100-200 MHz
• drive_bounds: Control bounds for [Ωᵣ(qubit), Ωᵢ(qubit), αᵣ(cavity), αᵢ(cavity)]
• multiply_by_2π: Whether to multiply by 2π (default false, assuming parameters already include it)

Example

# Standard cQED parameters
sys = TransmonCavitySystem(
    qubit_levels=4,
    cavity_levels=15,
    χ=2π * 32.8e-6,   # 32.8 kHz dispersive shift
    K_q=2π * 193e-3/2, # ~193 MHz anharmonicity
)

References

• Blais et al., "Circuit quantum electrodynamics," Rev. Mod. Phys. 93, 025005 (2021)
• Koch et al., "Charge-insensitive qubit design derived from Cooper pair box," Phys. Rev. A 76, 042319 (2007)

TransmonDipoleCoupling(
    g_ij::Float64,
    pair::Tuple{Int, Int},
    subsystem_levels::Vector{Int};
    lab_frame::Bool=false,
) -> QuantumSystemCoupling

TransmonDipoleCoupling(
    g_ij::Float64,
    pair::Tuple{Int, Int},
    sub_systems::Vector{QuantumSystem};
    kwargs...
) -> QuantumSystemCoupling

Returns a QuantumSystemCoupling object for a transmon qubit. In the lab frame, the Hamiltonian coupling term is

\[H = g_{ij} (a_i + a_i^\dagger) (a_j + a_j^\dagger)\]

In the rotating frame, the Hamiltonian coupling term is

\[H = g_{ij} (a_i a_j^\dagger + a_i^\dagger a_j)\]

where a_i is the annihilation operator for the ith transmon.

TransmonSystem(;
    ω::Float64=4.4153,  # GHz
    δ::Float64=0.17215, # GHz
    levels::Int=3,
    lab_frame::Bool=false,
    frame_ω::Float64=ω,
) -> QuantumSystem

Returns a QuantumSystem object for a transmon qubit, with the Hamiltonian

\[H = \omega a^\dagger a - \frac{\delta}{2} a^\dagger a^\dagger a a\]

where a is the annihilation operator.

Keyword Arguments

• ω: The frequency of the transmon, in GHz.
• δ: The anharmonicity of the transmon, in GHz.
• levels: The number of levels in the transmon.
• lab_frame: Whether to use the lab frame Hamiltonian, or an ω-rotating frame.
• frame_ω: The frequency of the rotating frame, in GHz.
• mutiply_by_2π: Whether to multiply the Hamiltonian by 2π, set to true by default because the frequency is in GHz.
• lab_frame_type: The type of lab frame Hamiltonian to use, one of (:duffing, :quartic, :cosine).
• drives: Whether to include drives in the Hamiltonian.

coherent_ket(α, levels)

Construct a coherent state $|α⟩$ in the Fock basis truncated to levels.

get_cat_controls(sys::AbstractQuantumSystem, α, N)

Compute static cat qubit controls for N time steps at coherent amplitude α, reading g2 and χ_aa from sys.global_params.

Returns a 2 × N matrix where row 1 is the buffer drive and row 2 is the Kerr correction drive. These are the steady-state controls that maintain a coherent state $|α⟩$ in the cat mode.

Arguments

• sys: A quantum system with g2 and χ_aa in its global_params
• α: Coherent state amplitude
• N: Number of time steps

lift(x::Char, i::Int, N::Int)::String

Embed a character into a string of the form 'I' * N at a specific position (meant for use with Piccolo.QuantumObjectUtils.operator_from_string).

is_linearly_dependent(M::AbstractMatrix; eps=eps(Float32), verbose=true)

Check if the columns of the matrix M are linearly dependent.

is_reachable(gate::AbstractMatrix{<:Number}, system::AbstractQuantumSystem; kwargs...)

Check if the gate is reachable using the given system.

Keyword Arguments

• use_drift::Bool=true: include drift Hamiltonian in the generators
• kwargs...: keyword arguments for is_reachable

is_reachable(gate, hamiltonians; kwargs...)

Check if the gate is reachable using the given hamiltonians.

Arguments

• gate::AbstractMatrix: target gate
• hamiltonians::AbstractVector{<:AbstractMatrix}: generators of the Lie algebra

Keyword Arguments

• subspace::AbstractVector{<:Int}=1:size(gate, 1): subspace indices
• compute_basis::Bool=true: compute the basis or use the Hamiltonians directly
• remove_trace::Bool=true: remove trace from generators
• verbose::Bool=true: print information about the operator algebra
• atol::Float32=eps(Float32): absolute tolerance

See also QuantumSystemUtils.operator_algebra.

operator_algebra(generators; kwargs...)

Compute the Lie algebra basis for the given generators.

Arguments

• generators::Vector{<:AbstractMatrix}: generators of the Lie algebra

Keyword Arguments

• return_layers::Bool=false: return the Lie tree layers
• normalize::Bool=false: normalize the basis
• verbose::Bool=false: print information
• remove_trace::Bool=true: remove trace from generators

A constant dictionary GATES containing common quantum gate matrices as complex-valued matrices. Each gate is represented by its unitary matrix.

• GATES[:I] - Identity: Leaves the state unchanged.
• GATES[:X] - Pauli-X (NOT): Flips the qubit state.
• GATES[:Y] - Pauli-Y: Rotates the qubit state around the Y-axis of the Bloch sphere.
• GATES[:Z] - Pauli-Z: Flips the phase of the qubit state.
• GATES[:H] - Hadamard: Creates superposition by transforming basis states.
• GATES[:CX] - Controlled-X (CNOT): Flips the 2nd qubit (target) if the first qubit (control) is |1⟩.
• GATES[:CZ] - Controlled-Z (CZ): Flips the phase of the 2nd qubit (target) if the 1st qubit (control) is |1⟩.
• GATES[:XI] - Complex: A gate for complex operations.
• GATES[:sqrtiSWAP] - Square root of iSWAP: Partially swaps two qubits with a phase.

The 2×2 Pauli matrics and identity.

annihilate(levels::Int)

Get the annihilation operator for a system with levels.

create(levels::Int)

Get the creation operator for a system with levels.

haar_identity(n::Int, radius::Number)

Generate a random unitary matrix close to the identity matrix using the Haar measure for an n-dimensional system with a given radius. The smaller the radius, the closer the matrix will be to the identity.

haar_random(n::Int)

Generate a random unitary matrix using the Haar measure for an n-dimensional system.

ket_from_bitstring(ket::String)

Get the state vector for a qubit system given a ket string ket of 0s and 1s.

ket_from_string(
    ket::String,
    levels::Vector{Int};
    level_dict=Dict(:g => 0, :e => 1, :f => 2, :h => 3, :i => 4, :j => 5, :k => 6, :l => 7),
    return_states=false
)

Construct a quantum state from a string ket representation.

operator_from_string(operator::String; lookup=PAULIS)

Reduce the string (each character is one key) via operators from a dictionary.

NamedTrajectory(qtraj::MultiKetTrajectory; kwargs...)
NamedTrajectory(qtraj::MultiKetTrajectory, N::Int; kwargs...)
NamedTrajectory(qtraj::MultiKetTrajectory, times::AbstractVector; kwargs...)

Convert an MultiKetTrajectory to a NamedTrajectory for optimization.

Stored Variables

• ψ̃1, ψ̃2, ...: Isomorphisms of each ket state
• u (or custom drive_name): Control values sampled at times
• du: Control derivatives (only for CubicSplinePulse)
• t: Times

Arguments

• N_or_times: One of:
  • nothing (default): Use native knot times from spline pulse
  • N::Int: Number of uniformly spaced time points
  • times::AbstractVector: Specific times to sample at

Keyword Arguments

• Δt_bounds: Optional tuple (lower, upper) for timestep bounds. If provided, enables free-time optimization (minimum-time problems). Default: nothing (no bounds).
• global_data: Optional Dict mapping global variable names to initial values (as vectors). Note: global variables are optimization variables without explicit box constraints.

NamedTrajectory(qtraj::DensityTrajectory; kwargs...)
NamedTrajectory(qtraj::DensityTrajectory, N::Int; kwargs...)
NamedTrajectory(qtraj::DensityTrajectory, times::AbstractVector; kwargs...)

Convert a DensityTrajectory to a NamedTrajectory for optimization.

Uses the compact density isomorphism (n² real parameters) which exploits the Hermiticity of density matrices. This halves the state dimension compared to the full [vec(Re(ρ)); vec(Im(ρ))] representation.

Stored Variables

• ρ⃗̃: Compact isomorphism of the density matrix (n² real parameters)
• u (or custom drive_name): Control values sampled at times
• du: Control derivatives (only for CubicSplinePulse)
• t: Times

Arguments

• N_or_times: One of:
  • nothing (default): Use native knot times from spline pulse
  • N::Int: Number of uniformly spaced time points
  • times::AbstractVector: Specific times to sample at

Keyword Arguments

• Δt_bounds: Optional tuple (lower, upper) for timestep bounds. If provided, enables free-time optimization (minimum-time problems). Default: nothing (no bounds).
• global_data: Optional Dict mapping global variable names to initial values (as vectors). Note: global variables are optimization variables without explicit box constraints.

NamedTrajectory(qtraj::KetTrajectory; kwargs...)
NamedTrajectory(qtraj::KetTrajectory, N::Int; kwargs...)
NamedTrajectory(qtraj::KetTrajectory, times::AbstractVector; kwargs...)

Convert a KetTrajectory to a NamedTrajectory for optimization.

Stored Variables

• ψ̃: Isomorphism of ket state (real representation)
• u (or custom drive_name): Control values sampled at times
• du: Control derivatives (only for CubicSplinePulse)
• t: Times

Arguments

• N_or_times: One of:
  • nothing (default): Use native knot times from spline pulse
  • N::Int: Number of uniformly spaced time points
  • times::AbstractVector: Specific times to sample at

Keyword Arguments

• Δt_bounds: Optional tuple (lower, upper) for timestep bounds. If provided, enables free-time optimization (minimum-time problems). Default: nothing (no bounds).
• global_data: Optional Dict mapping global variable names to initial values (as vectors). Note: global variables are optimization variables without explicit box constraints.

NamedTrajectory(qtraj::UnitaryTrajectory; kwargs...)
NamedTrajectory(qtraj::UnitaryTrajectory, N::Int; kwargs...)
NamedTrajectory(qtraj::UnitaryTrajectory, times::AbstractVector; kwargs...)

Convert a UnitaryTrajectory to a NamedTrajectory for optimization.

The trajectory stores actual times :t (not timesteps :Δt), which is required for time-dependent integrators used with SplinePulseProblem.

Stored Variables

• Ũ⃗: Isomorphism of unitary (vectorized real representation)
• u (or custom drive_name): Control values sampled at times
• du: Control derivatives (only for CubicSplinePulse)
• t: Times

Arguments

• qtraj: The quantum trajectory to convert
• N_or_times: One of:
  • nothing (default): Use native knot times from spline pulse (error for non-spline pulses)
  • N::Int: Number of uniformly spaced time points
  • times::AbstractVector: Specific times to sample at

Keyword Arguments

• Δt_bounds: Optional tuple (lower, upper) for timestep bounds. If provided, enables free-time optimization (minimum-time problems). Default: nothing (no bounds).
• global_data: Optional Dict mapping global variable names to initial values (as vectors). Note: global variables are optimization variables without explicit box constraints.

Returns

A NamedTrajectory suitable for direct collocation optimization.

NamedTrajectory(sampling::SamplingTrajectory, N::Int)
NamedTrajectory(sampling::SamplingTrajectory, times::AbstractVector)

Convert a SamplingTrajectory to a NamedTrajectory for optimization.

Creates a trajectory with multiple state variables (one per system), all sharing the same control pulse. Each state gets a numeric suffix:

• UnitaryTrajectory base → :Ũ⃗1, :Ũ⃗2, ...
• KetTrajectory base → :ψ̃1, :ψ̃2, ...

For robust optimization, each state variable represents the evolution under a different system (e.g., parameter variations), but all share the same controls.

Example

# Create sampling trajectory with 3 system variations
sampling = SamplingTrajectory(base_qtraj, [sys1, sys2, sys3])

# Convert to NamedTrajectory with 51 timesteps
traj = NamedTrajectory(sampling, 51)
# Result has: :Ũ⃗1, :Ũ⃗2, :Ũ⃗3, :u, :Δt, :t

Keyword Arguments

• Δt_bounds: Optional tuple (lower, upper) for timestep bounds. If provided, enables free-time optimization (minimum-time problems). Default: nothing (no bounds).

AbstractQuantumTrajectory{P<:AbstractPulse}

Abstract type for quantum trajectories that wrap physics (system, pulse, solution, goal). Parametric on pulse type P to enable dispatch in problem templates.

All concrete subtypes should implement:

• state_name(traj) - Get the state variable symbol (fixed per type)
• drive_name(traj) - Get the drive variable symbol (from pulse)
• time_name(traj) - Get the time variable symbol (fixed :t)
• timestep_name(traj) - Get the timestep variable symbol (fixed :Δt)
• duration(traj) - Get the duration (from pulse)

DensityTrajectory{P<:AbstractPulse, S<:ODESolution} <: AbstractQuantumTrajectory{P}

Trajectory for open quantum systems (Lindblad dynamics).

Fields

• system::OpenQuantumSystem: The open quantum system
• pulse::P: The control pulse
• initial::Matrix{ComplexF64}: Initial density matrix ρ₀
• goal::Matrix{ComplexF64}: Target density matrix ρ_goal
• solution::S: Pre-computed ODE solution

Callable

traj(t) returns the density matrix at time t.

DensityTrajectory(system, initial, goal, T::Real; drive_name=:u, algorithm=Tsit5(), abstol=1e-8, reltol=1e-8)

Convenience constructor that creates a zero pulse of duration T.

DensityTrajectory(system, pulse, initial, goal; algorithm=Tsit5(), abstol=1e-8, reltol=1e-8, n_save=101)

Create a density matrix trajectory by solving the Lindblad master equation.

Arguments

• system::OpenQuantumSystem: The open quantum system
• pulse::AbstractPulse: The control pulse
• initial::Matrix: Initial density matrix ρ₀
• goal::Matrix: Target density matrix ρ_goal

Keyword Arguments

• algorithm: ODE solver algorithm (default: Tsit5())
• abstol: Absolute tolerance (default: 1e-8)
• reltol: Relative tolerance (default: 1e-8)
• n_save: Number of output time points (default: 101)

KetTrajectory{P<:AbstractPulse, S<:ODESolution} <: AbstractQuantumTrajectory{P}

Trajectory for quantum state transfer. The ODE solution is computed at construction.

Fields

• system::QuantumSystem: The quantum system
• pulse::P: The control pulse
• initial::Vector{ComplexF64}: Initial state |ψ₀⟩
• goal::Vector{ComplexF64}: Target state |ψ_goal⟩
• solution::S: Pre-computed ODE solution

Callable

traj(t) returns the state at time t by interpolating the solution.

KetTrajectory(system, pulse, initial, goal; algorithm=MagnusAdapt4(), abstol=1e-8, reltol=1e-8, n_save=101)

Create a ket trajectory by solving the Schrödinger equation.

Arguments

• system::QuantumSystem: The quantum system
• pulse::AbstractPulse: The control pulse
• initial::Vector: Initial state |ψ₀⟩
• goal::Vector: Target state |ψ_goal⟩

Keyword Arguments

• algorithm: ODE solver algorithm (default: MagnusAdapt4())
• abstol: Absolute tolerance for adaptive integration (default: 1e-8)
• reltol: Relative tolerance for adaptive integration (default: 1e-8)
• n_save: Number of output time points (default: 101)

KetTrajectory(system, initial, goal, T::Real; drive_name=:u, algorithm=MagnusAdapt4(), abstol=1e-8, reltol=1e-8)

Convenience constructor that creates a zero pulse of duration T.

Arguments

• system::QuantumSystem: The quantum system
• initial::Vector: Initial state |ψ₀⟩
• goal::Vector: Target state |ψ_goal⟩
• T::Real: Duration of the pulse

Keyword Arguments

• drive_name::Symbol: Name of the drive variable (default: :u)
• algorithm: ODE solver algorithm (default: MagnusAdapt4())
• abstol: Absolute tolerance (default: 1e-8)
• reltol: Relative tolerance (default: 1e-8)

MultiKetTrajectory{P<:AbstractPulse, S} <: AbstractQuantumTrajectory{P}

Trajectory for multi-state transfer with a shared pulse. Useful for state-to-state problems with multiple initial/goal pairs.

Fields

• system::QuantumSystem: The quantum system
• pulse::P: The shared control pulse
• initials::Vector{Vector{ComplexF64}}: Initial states
• goals::Vector{Vector{ComplexF64}}: Target states
• weights::Vector{Float64}: Weights for fidelity calculation
• solution::S: Pre-computed ensemble solution

Callable

traj(t) returns a vector of states at time t. traj[i] returns the i-th trajectory's solution.

MultiKetTrajectory(system, pulse, initials, goals; weights=..., algorithm=MagnusAdapt4(), abstol=1e-8, reltol=1e-8, n_save=101)

Create a multi-ket trajectory by solving multiple Schrödinger equations.

Arguments

• system::QuantumSystem: The quantum system
• pulse::AbstractPulse: The shared control pulse
• initials::Vector{Vector}: Initial states
• goals::Vector{Vector}: Target states

Keyword Arguments

• weights: Weights for fidelity (default: uniform)
• algorithm: ODE solver algorithm (default: MagnusAdapt4())
• abstol: Absolute tolerance for adaptive integration (default: 1e-8)
• reltol: Relative tolerance for adaptive integration (default: 1e-8)
• n_save: Number of output time points (default: 101)

MultiKetTrajectory(system, initials, goals, T::Real; weights=..., drive_name=:u, algorithm=MagnusAdapt4(), abstol=1e-8, reltol=1e-8)

Convenience constructor that creates a zero pulse of duration T.

Arguments

• system::QuantumSystem: The quantum system
• initials::Vector{Vector}: Initial states
• goals::Vector{Vector}: Target states
• T::Real: Duration of the pulse

Keyword Arguments

• weights: Weights for fidelity (default: uniform)
• drive_name::Symbol: Name of the drive variable (default: :u)
• algorithm: ODE solver algorithm (default: MagnusAdapt4())
• abstol: Absolute tolerance (default: 1e-8)
• reltol: Relative tolerance (default: 1e-8)

SamplingTrajectory{QT<:AbstractQuantumTrajectory} <: AbstractQuantumTrajectory

Wrapper for robust optimization over multiple systems with shared controls.

Used for sampling-based robust optimization where:

• All systems share the same control pulse
• Each system has different dynamics (e.g., parameter variations)
• Optimization minimizes weighted fidelity across all systems

This type does NOT store a NamedTrajectory - use NamedTrajectory(sampling, N) for conversion.

Fields

• base_trajectory::QT: Base quantum trajectory (defines pulse, initial, goal)
• systems::Vector{<:AbstractQuantumSystem}: Multiple systems to optimize over
• weights::Vector{Float64}: Weights for each system in objective

Example

sys_nom = QuantumSystem(...)
sys_variations = [QuantumSystem(...) for _ in 1:3]  # Parameter variations
qtraj = UnitaryTrajectory(sys_nom, pulse, U_goal)
sampling = SamplingTrajectory(qtraj, sys_variations, [0.5, 0.3, 0.2])

# Convert to NamedTrajectory for optimization
traj = NamedTrajectory(sampling, 51)  # Creates :Ũ⃗1, :Ũ⃗2, :Ũ⃗3

SamplingTrajectory(base_trajectory, systems; weights=nothing)

Create a SamplingTrajectory for robust optimization.

Arguments

• base_trajectory: Base quantum trajectory (defines pulse, initial, goal)
• systems: Vector of systems with parameter variations

Keyword Arguments

• weights: Optional weights for each system (default: equal weights)

UnitaryTrajectory{P<:AbstractPulse, S<:ODESolution, G} <: AbstractQuantumTrajectory{P}

Trajectory for unitary gate synthesis. The ODE solution is computed at construction.

Fields

• system::QuantumSystem: The quantum system
• pulse::P: The control pulse (stores drive_name)
• initial::Matrix{ComplexF64}: Initial unitary (default: identity)
• goal::G: Target unitary operator (AbstractPiccoloOperator or Matrix)
• solution::S: Pre-computed ODE solution

Callable

traj(t) returns the unitary at time t by interpolating the solution.

Conversion to NamedTrajectory

Use NamedTrajectory(traj, N) or NamedTrajectory(traj, times) for optimization.

UnitaryTrajectory(system, pulse, goal; initial=I, algorithm=MagnusAdapt4(), abstol=1e-8, reltol=1e-8, n_save=101)

Create a unitary trajectory by solving the Schrödinger equation.

Arguments

• system::QuantumSystem: The quantum system
• pulse::AbstractPulse: The control pulse
• goal: Target unitary (Matrix or AbstractPiccoloOperator)

Keyword Arguments

• initial: Initial unitary (default: identity matrix)
• algorithm: ODE solver algorithm (default: MagnusAdapt4())
• abstol: Absolute tolerance for adaptive integration (default: 1e-8)
• reltol: Relative tolerance for adaptive integration (default: 1e-8)
• n_save: Number of output time points for plotting/interpolation (default: 101)

UnitaryTrajectory(system, goal, T::Real; drive_name=:u, algorithm=MagnusAdapt4(), abstol=1e-8, reltol=1e-8)

Convenience constructor that creates a zero pulse of duration T.

Arguments

• system::QuantumSystem: The quantum system
• goal: Target unitary (Matrix or AbstractPiccoloOperator)
• T::Real: Duration of the pulse

Keyword Arguments

• drive_name::Symbol: Name of the drive variable (default: :u)
• algorithm: ODE solver algorithm (default: MagnusAdapt4())
• abstol: Absolute tolerance (default: 1e-8)
• reltol: Relative tolerance (default: 1e-8)

drive_name(traj::AbstractQuantumTrajectory)

Get the drive/control variable name from the trajectory's pulse.

duration(traj)

Get the duration of a trajectory (from its pulse).

_add_global_data_to_kwargs(nt_kwargs, global_data)

Helper function to process global variables and add them to NamedTrajectory kwargs. Converts Dict{Symbol, Vector} to flat vector and components NamedTuple.

Arguments

• nt_kwargs: Existing NamedTuple of kwargs to merge with
• global_data: Dict mapping global variable names to vectors of values

Returns

Merged NamedTuple with globaldata and globalcomponents added

_get_control_data(pulse::CubicSplinePulse, times, sys)

For CubicSplinePulse: return u and du data with system bounds and boundary conditions. Uses the pulse's drive_name to determine variable naming.

When times matches the pulse's native knot times, extracts stored u and du directly. When resampling to different times, samples u via interpolation and computes du via ForwardDiff to get the true spline derivative.

Returns

• data: NamedTuple with control data
• control_names: Tuple of control variable names
• bounds: NamedTuple with control bounds
• initial_constraints: NamedTuple with initial value constraints
• final_constraints: NamedTuple with final value constraints

_get_control_data(pulse::GaussianPulse, times, sys)

For GaussianPulse: sample as u values with system bounds and boundary conditions. Uses the pulse's drive_name to determine variable naming.

Returns

• data: NamedTuple with control data
• control_names: Tuple of control variable names
• bounds: NamedTuple with control bounds
• initial_constraints: NamedTuple with initial value constraints (empty for parametric pulses)
• final_constraints: NamedTuple with final value constraints (empty for parametric pulses)

_get_control_data(pulse::Union{ZeroOrderPulse, LinearSplinePulse}, times, sys)

For ZeroOrderPulse and LinearSplinePulse: return u data with system bounds and boundary conditions. Uses the pulse's drive_name to determine variable naming.

Returns

• data: NamedTuple with control data
• control_names: Tuple of control variable names
• bounds: NamedTuple with control bounds
• initial_constraints: NamedTuple with initial value constraints
• final_constraints: NamedTuple with final value constraints

_get_drive_bounds(sys::QuantumSystem)

Extract drive bounds from system as tuple of (lower, upper) vectors.

_named_tuple(pairs...)

Create a NamedTuple from pairs of (Symbol, value). This is needed when keys are dynamic (stored in variables).

Example: name = :x namedtuple(name => 1, :y => 2) # Returns (x = 1, y = 2)

_sample_times(traj, times::AbstractVector)

Return times as a Float64 vector.

_sample_times(traj, N::Int)

Generate N uniformly spaced times for sampling.

_sample_times(traj, ::Nothing)

For spline pulses, extract native knot times. For other pulses, error.

extract_pulse(qtraj::AbstractQuantumTrajectory, traj::NamedTrajectory)

Extract an optimized pulse from a NamedTrajectory.

This function extracts the control values from the optimized trajectory and creates a new pulse object of the same type as the original pulse in qtraj.

The extraction process depends on the pulse type:

• ZeroOrderPulse, LinearSplinePulse: Extracts u (drive variable)
• CubicSplinePulse: Extracts both u and du (derivative variable)

Arguments

• qtraj: Original quantum trajectory (provides pulse type and drive names)
• traj: Optimized NamedTrajectory with new control values

Returns

A new pulse of the same type as qtraj.pulse with optimized control values.

Example

# After optimization
solve!(prob)
new_pulse = extract_pulse(qtraj, prob.trajectory)
rollout!(qtraj, new_pulse)

get_goal(traj)

Get the goal state/operator from a trajectory.

get_initial(traj)

Get the initial state/operator from a trajectory.

get_pulse(traj)

Get the control pulse from a trajectory.

get_solution(traj)

Get the ODE solution from a trajectory.

get_system(traj)

Get the quantum system from a trajectory.

get_systems(sampling::SamplingTrajectory)

Get all systems in the sampling trajectory.

get_weights(sampling::SamplingTrajectory)

Get the weights for each system.

state_name(::AbstractQuantumTrajectory)

Get the fixed state variable name for a trajectory type.

• UnitaryTrajectory → :Ũ⃗
• KetTrajectory → :ψ̃
• MultiKetTrajectory → :ψ̃ (with index appended: :ψ̃1, :ψ̃2, etc.)
• DensityTrajectory → :ρ⃗̃

state_names(traj::MultiKetTrajectory)

Get all state names for an ensemble trajectory (:ψ̃1, :ψ̃2, etc.)

state_names(sampling::SamplingTrajectory)

Get the state variable names for all systems (e.g., [:Ũ⃗1, :Ũ⃗2, :Ũ⃗3]).

time_name(::AbstractQuantumTrajectory)

Get the time variable name (always :t).

timestep_name(::AbstractQuantumTrajectory)

Get the timestep variable name (always :Δt).

Rollouts._update_system!(qtraj::DensityTrajectory, sys::OpenQuantumSystem)

Update the system field in a DensityTrajectory with a new OpenQuantumSystem (typically with updated global parameters after optimization).

Rollouts._update_system!(qtraj::KetTrajectory, sys::QuantumSystem)

Update the system field in a KetTrajectory with a new QuantumSystem (typically with updated global parameters after optimization).

Rollouts._update_system!(qtraj::MultiKetTrajectory, sys::QuantumSystem)

Update the system field in a MultiKetTrajectory with a new QuantumSystem (typically with updated global parameters after optimization).

Rollouts._update_system!(qtraj::SamplingTrajectory, sys::QuantumSystem)

Update the system in the base_trajectory of a SamplingTrajectory. Note: This only updates the base trajectory's system, not the systems array. For updating parameter variations in the systems array, that should be done through the SamplingTrajectory constructor or direct modification.

Rollouts._update_system!(qtraj::UnitaryTrajectory, sys::QuantumSystem)

Update the system field in a UnitaryTrajectory with a new QuantumSystem (typically with updated global parameters after optimization).

fidelity(traj::DensityTrajectory)

Compute the fidelity between the final density matrix and the goal. Uses trace fidelity: F = tr(ρfinal * ρgoal)

fidelity(traj::KetTrajectory)

Compute the fidelity between the final state and the goal.

fidelity(traj::MultiKetTrajectory)

Compute the weighted average fidelity across all state transfers.

fidelity(traj::SamplingTrajectory; kwargs...)

Compute the fidelity for each system in the sampling trajectory.

Returns a vector of fidelities, one per system, by rolling out the current pulse with each system and computing the fidelity against the goal.

fidelity(traj::UnitaryTrajectory; subspace=nothing, phases=nothing)

Compute the fidelity between the final unitary and the goal.

For EmbeddedOperator goals, uses the Pedersen average gate fidelity on the computational subspace, matching the metric the optimizer minimizes:

F = 1/(n(n+1)) * (|tr(M'M)| + |tr(M)|²),   M = U_goal' * U_sub

For standard goals (or when an explicit subspace is given), uses the standard unitary fidelity |tr(U'U_goal)|²/N².

When phases is provided (a vector of per-qubit Z-phases of length n_qubits), the goal on the computational subspace is adjusted by Diagonal(phase_vec) * U_goal_sub before computing fidelity, where phase_vec is a length-n_sub vector whose entries are products of exp(im * θⱼ) for each qubit j that is in the |1⟩ state of the corresponding computational basis vector (via binary decomposition of the basis index). This matches the free-phase objective used during optimization.

rollout!(qtraj::DensityTrajectory, pulse::AbstractPulse; algorithm=Tsit5(), n_save=101, abstol=1e-8, reltol=1e-8)

Update density trajectory in-place with a new pulse. Note: Default algorithm is Tsit5() since density evolution uses standard ODE solvers. See rollout!(::UnitaryTrajectory, ::AbstractPulse) for details.

rollout!(qtraj::DensityTrajectory; algorithm=Tsit5(), n_save=101, abstol=1e-8, reltol=1e-8, kwargs...)

Update density trajectory in-place with same pulse but different ODE parameters. Note: Default algorithm is Tsit5() since density evolution uses standard ODE solvers. See rollout!(::UnitaryTrajectory; kwargs...) for details.

rollout!(qtraj::KetTrajectory, pulse::AbstractPulse; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8)

Update ket trajectory in-place with a new pulse. See rollout!(::UnitaryTrajectory, ::AbstractPulse) for details.

rollout!(qtraj::KetTrajectory; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8, kwargs...)

Update ket trajectory in-place with same pulse but different ODE parameters. See rollout!(::UnitaryTrajectory; kwargs...) for details.

rollout!(qtraj::MultiKetTrajectory, pulse::AbstractPulse; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8)

Update multi-ket trajectory in-place with a new pulse. See rollout!(::UnitaryTrajectory, ::AbstractPulse) for details.

rollout!(qtraj::MultiKetTrajectory; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8, kwargs...)

Update multi-ket trajectory in-place with same pulse but different ODE parameters. See rollout!(::UnitaryTrajectory; kwargs...) for details.

rollout!(qtraj::SamplingTrajectory, pulse::AbstractPulse; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8)

Update sampling trajectory's base trajectory in-place with a new pulse. Delegates to the base trajectory's rollout! method.

rollout!(qtraj::SamplingTrajectory; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8, kwargs...)

Update sampling trajectory's base trajectory in-place with new ODE parameters. Delegates to the base trajectory's rollout! method.

rollout!(qtraj::UnitaryTrajectory, pulse::AbstractPulse; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8)

Update quantum trajectory in-place with a new pulse by re-solving the ODE. Mutates qtraj.pulse and qtraj.solution.

Arguments

• qtraj::UnitaryTrajectory: The trajectory to update
• pulse::AbstractPulse: The new control pulse

Keyword Arguments

• algorithm: ODE solver algorithm (default: MagnusAdapt4())
• n_save::Int: Number of output time points (default: 101)
• abstol: Absolute tolerance (default: 1e-8)
• reltol: Relative tolerance (default: 1e-8)

Example

qtraj = UnitaryTrajectory(sys, old_pulse, goal)
rollout!(qtraj, new_pulse)  # Updates qtraj in-place
fid = fidelity(qtraj)  # Uses new solution

See also: rollout

rollout!(qtraj::UnitaryTrajectory; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8, kwargs...)

Update quantum trajectory in-place by re-solving with same pulse but different ODE parameters. Mutates qtraj.solution.

Useful for comparing different solvers or tolerances.

Keyword Arguments

• algorithm: ODE solver algorithm (default: MagnusAdapt4())
• n_save::Int: Number of output time points (default: 101)
• abstol: Absolute tolerance (default: 1e-8)
• reltol: Relative tolerance (default: 1e-8)
• Additional kwargs passed to solve

Example

qtraj = UnitaryTrajectory(sys, pulse, goal)

# Compare Magnus vs Runge-Kutta
rollout!(qtraj; algorithm=Tsit5(), abstol=1e-10)
fid_rk = fidelity(qtraj)

See also: rollout

rollout(qtraj::DensityTrajectory, pulse::AbstractPulse; algorithm=Tsit5(), n_save=101, abstol=1e-8, reltol=1e-8)

Create a new density trajectory by rolling out a new pulse. Note: Default algorithm is Tsit5() since density evolution uses standard ODE solvers. See rollout(::UnitaryTrajectory, ::AbstractPulse) for details.

rollout(qtraj::DensityTrajectory; algorithm=Tsit5(), n_save=101, abstol=1e-8, reltol=1e-8, kwargs...)

Re-solve density trajectory with same pulse but different ODE parameters. Note: Default algorithm is Tsit5() since density evolution uses standard ODE solvers. See rollout(::UnitaryTrajectory; kwargs...) for details.

rollout(qtraj::KetTrajectory, pulse::AbstractPulse; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8)

Create a new ket trajectory by rolling out a new pulse. See rollout(::UnitaryTrajectory, ::AbstractPulse) for details.

rollout(qtraj::KetTrajectory; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8, kwargs...)

Re-solve ket trajectory with same pulse but different ODE parameters. See rollout(::UnitaryTrajectory; kwargs...) for details.

rollout(qtraj::MultiKetTrajectory, pulse::AbstractPulse; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8)

Create a new multi-ket trajectory by rolling out a new pulse. See rollout(::UnitaryTrajectory, ::AbstractPulse) for details.

rollout(qtraj::MultiKetTrajectory; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8, kwargs...)

Re-solve multi-ket trajectory with same pulse but different ODE parameters. See rollout(::UnitaryTrajectory; kwargs...) for details.

rollout(qtraj::UnitaryTrajectory, pulse::AbstractPulse; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8)

Create a new quantum trajectory by rolling out a new pulse through the system. Returns a new UnitaryTrajectory with the updated pulse and solution.

Arguments

• qtraj::UnitaryTrajectory: The base trajectory (provides system, initial, goal)
• pulse::AbstractPulse: The new control pulse to roll out

Keyword Arguments

• algorithm: ODE solver algorithm (default: MagnusAdapt4())
• n_save::Int: Number of output time points (default: 101)
• abstol: Absolute tolerance (default: 1e-8)
• reltol: Relative tolerance (default: 1e-8)

Example

qtraj = UnitaryTrajectory(sys, old_pulse, goal)

# Roll out a new pulse
qtraj_new = rollout(qtraj, new_pulse)

# Check fidelity
fid = fidelity(qtraj_new)

See also: extract_pulse, rollout!, fidelity

rollout(qtraj::UnitaryTrajectory; algorithm=MagnusAdapt4(), n_save=101, abstol=1e-8, reltol=1e-8, kwargs...)

Re-solve the trajectory with the same pulse but different ODE parameters. Returns a new UnitaryTrajectory with the updated solution.

Useful for comparing different solvers or tolerances.

Keyword Arguments

• algorithm: ODE solver algorithm (default: MagnusAdapt4())
• n_save::Int: Number of output time points (default: 101)
• abstol: Absolute tolerance (default: 1e-8)
• reltol: Relative tolerance (default: 1e-8)
• Additional kwargs passed to solve

Example

qtraj = UnitaryTrajectory(sys, pulse, goal)

# Compare Magnus vs Runge-Kutta
qtraj_rk = rollout(qtraj; algorithm=Tsit5(), abstol=1e-10)
fid_magnus = fidelity(qtraj)
fid_rk = fidelity(qtraj_rk)

See also: rollout!

AbstractPulse

Abstract type for all pulse types. All pulses are callable: pulse(t) returns the control vector at time t.

AbstractSplinePulse <: AbstractPulse

Abstract type for spline-based pulses (linear and cubic interpolation). These pulses use the spline coefficients as optimization variables.

CompositePulse(pulses::Vector{<:AbstractPulse}, mode::Symbol=:interleave)

Create a composite pulse from multiple component pulses.

Arguments

• pulses: Vector of pulse objects to combine
• mode: How to combine the drives
  • :interleave - Interleave drives: [p1d1, p2d1, p1d2, p2d2, ...]
  • :concatenate - Concatenate drives: [p1d1, p1d2, ..., p2d1, p2d2, ...]

Example

# For MS gate with 2 ions: [Ω₁, φ₁, Ω₂, φ₂]
Ω_pulse = GaussianPulse([Ω₁, Ω₂], σ, T)  # 2 drives
φ_pulse = ErfPulse([φ₁, φ₂], σ, T)        # 2 drives
pulse = CompositePulse([Ω_pulse, φ_pulse], :interleave)
# Result: pulse(t) = [Ω₁(t), φ₁(t), Ω₂(t), φ₂(t)]

CompositePulse{F<:Function} <: AbstractPulse

Composite pulse that combines multiple pulse objects by interleaving their drives.

Useful for creating pulses with different shapes for different control types, such as Gaussian amplitude + erf phase for trapped ion gates.

Fields

• f::F: Function that evaluates the composite pulse
• pulses::Vector{<:AbstractPulse}: Component pulses
• drive_mapping::Vector{Vector{Int}}: Maps pulse i, drive j to composite drive index
• duration::Float64: Total pulse duration (must match for all components)
• n_drives::Int: Total number of drives across all pulses

Example

# Amplitude: Gaussian (2 drives for 2 ions)
Ω_pulse = GaussianPulse([Ω_max, Ω_max], σ, T)

# Phase: Error function (2 drives for 2 ions)
φ_pulse = ErfPulse([φ_max, φ_max], σ, T)

# Composite: [Ω₁, φ₁, Ω₂, φ₂] - interleaved
pulse = CompositePulse([Ω_pulse, φ_pulse], :interleave)

CubicSplinePulse{I<:CubicHermiteSpline} <: AbstractPulse

Pulse with cubic Hermite spline interpolation. Uses both control values AND derivatives for exact reconstruction after optimization.

Fields

• controls::I: CubicHermiteSpline from DataInterpolations
• duration::Float64: Total pulse duration
• n_drives::Int: Number of control drives
• drive_name::Symbol: Name of the drive variable (default :u)
• initial_value::Vector{Float64}: Initial boundary condition (default: zeros)
• final_value::Vector{Float64}: Final boundary condition (default: zeros)

CubicSplinePulse(controls::AbstractMatrix, derivatives::AbstractMatrix, times::AbstractVector; drive_name=:u, initial_value=nothing, final_value=nothing)

Create a cubic Hermite spline pulse from control values, derivatives, and times.

Arguments

• controls: Matrix of size (n_drives, n_times) with control values
• derivatives: Matrix of size (n_drives, n_times) with control derivatives
• times: Vector of sample times (must start at 0)

Keyword Arguments

• drive_name: Name of the drive variable (default :u)
• initial_value: Initial boundary condition (default: zeros(n_drives))
• final_value: Final boundary condition (default: zeros(n_drives))

CubicSplinePulse(controls::AbstractMatrix, times::AbstractVector; drive_name=:u, initial_value=nothing, final_value=nothing)

Create a cubic Hermite spline pulse with zero derivatives at all knot points. Useful for initial guesses where smoothness constraints will be enforced by optimizer.

Arguments

• controls: Matrix of size (n_drives, n_times) with control values
• times: Vector of sample times (must start at 0)

Keyword Arguments

• drive_name: Name of the drive variable (default :u)
• initial_value: Initial boundary condition (default: zeros(n_drives))
• final_value: Final boundary condition (default: zeros(n_drives))

CubicSplinePulse(pulse::AbstractPulse, n_samples::Int; kwargs...)
CubicSplinePulse(pulse::AbstractPulse, times::AbstractVector; kwargs...)

Convert any pulse to a CubicSplinePulse by sampling at specified times. Derivatives are computed using ForwardDiff for automatic differentiation.

Useful for initializing optimization problems with smooth analytic pulse shapes.

Arguments

• pulse: Source pulse (GaussianPulse, ErfPulse, CompositePulse, etc.)
• n_samples: Number of uniformly spaced samples (alternative to times)
• times: Specific sample times (alternative to n_samples)

Keyword Arguments

• drive_name: Name for the drive variable (default: :du)
• initial_value: Initial boundary condition (default: pulse(0.0))
• final_value: Final boundary condition (default: pulse(duration))

Example

gaussian = GaussianPulse([1.0, 2.0], 0.1, 1.0)
cubic = CubicSplinePulse(gaussian, 50)  # 50 samples with ForwardDiff derivatives

CubicSplinePulse(traj::NamedTrajectory; drive_name=:u, derivative_name=:du)

Construct a CubicSplinePulse (Hermite) from a NamedTrajectory using both control values and derivatives.

Arguments

• traj: NamedTrajectory with control and derivative data

Keyword Arguments

• drive_name: Name of the drive component (default: :u)
• derivative_name: Name of the derivative component (default: :du)

ErfPulse{F<:Function} <: AbstractPulse

Analytic error function pulse for phase compensation in trapped ion gates.

The error function profile is commonly used to compensate AC Stark shifts in Mølmer-Sørensen gates, where φ(t) ∝ erf(√2 (t - t₀)/σ) cancels time-varying phases from off-resonant spectator modes.

u_i(t) = amplitudes[i] * erf(√2 * (t - centers[i]) / sigmas[i])

Typically scaled to range [0, 1] or [-1, 1] by adjusting amplitude.

Fields

• f::F: Function that evaluates the pulse
• amplitudes::Vector{Float64}: Peak amplitude for each drive
• sigmas::Vector{Float64}: Width parameter for each drive
• centers::Vector{Float64}: Center time for each drive
• duration::Float64: Total pulse duration
• n_drives::Int: Number of control drives

References

• Mizrahi et al., "Realization and Calibration of Continuously Parameterized Two-Qubit Gates...", IEEE TQE (2024), Figure 7b

ErfPulse(amplitudes, sigmas, centers, duration)

Create an error function pulse with per-drive parameters.

Arguments

• amplitudes: Peak amplitude for each drive
• sigmas: Width parameter for each drive (controls steepness)
• centers: Center time for each drive (inflection point)
• duration: Total pulse duration

Example

using SpecialFunctions: erf

# Phase compensation for MS gate
φ_max = π/4  # Maximum phase shift
T = 50.0     # Gate duration
σ = T/4      # Width parameter

pulse = ErfPulse([φ_max], [σ], [T/2], T)

ErfPulse(amplitudes, sigma, duration; center=duration/2)

Create an error function pulse with shared sigma and center across all drives.

Arguments

• amplitudes: Peak amplitude for each drive
• sigma: Shared width parameter for all drives
• duration: Total pulse duration

Keyword Arguments

• center: Shared center time (default: duration/2)

GaussianPulse{F<:Function} <: AbstractPulse

Analytic Gaussian pulse. Each drive has its own amplitude, width (sigma), and center.

u_i(t) = amplitudes[i] * exp(-(t - centers[i])² / (2 * sigmas[i]²))

Fields

• f::F: Function that evaluates the pulse
• amplitudes::Vector{Float64}: Peak amplitude for each drive
• sigmas::Vector{Float64}: Gaussian width for each drive
• centers::Vector{Float64}: Center time for each drive
• duration::Float64: Total pulse duration
• n_drives::Int: Number of control drives

GaussianPulse(amplitudes, sigmas, centers, duration)

Create a Gaussian pulse with per-drive parameters.

Arguments

• amplitudes: Peak amplitude for each drive
• sigmas: Gaussian width (standard deviation) for each drive
• centers: Center time for each drive
• duration: Total pulse duration

GaussianPulse(amplitudes, sigma, duration; center=duration/2)

Create a Gaussian pulse with shared sigma and center across all drives.

Arguments

• amplitudes: Peak amplitude for each drive
• sigma: Shared Gaussian width for all drives
• duration: Total pulse duration

Keyword Arguments

• center: Shared center time (default: duration/2)

LinearSplinePulse{I<:LinearInterpolation} <: AbstractSplinePulse

Pulse with linear interpolation between sample points.

Fields

• controls::I: LinearInterpolation from DataInterpolations
• duration::Float64: Total pulse duration
• n_drives::Int: Number of control drives
• drive_name::Symbol: Name of the drive variable (default :u)
• initial_value::Vector{Float64}: Initial boundary condition (default: zeros)
• final_value::Vector{Float64}: Final boundary condition (default: zeros)

LinearSplinePulse(controls::AbstractMatrix, times::AbstractVector; drive_name=:u, initial_value=nothing, final_value=nothing)

Create a linearly interpolated pulse from control samples and times.

Arguments

• controls: Matrix of size (n_drives, n_times) with control values
• times: Vector of sample times (must start at 0)

Keyword Arguments

• drive_name: Name of the drive variable (default :u)
• initial_value: Initial boundary condition (default: zeros(n_drives))
• final_value: Final boundary condition (default: zeros(n_drives))

LinearSplinePulse(pulse::AbstractPulse, n_samples::Int; kwargs...)
LinearSplinePulse(pulse::AbstractPulse, times::AbstractVector; kwargs...)

Convert any pulse to a LinearSplinePulse by sampling at specified times.

Useful for initializing optimization problems with analytic pulse shapes.

Arguments

• pulse: Source pulse (GaussianPulse, ErfPulse, CompositePulse, etc.)
• n_samples: Number of uniformly spaced samples (alternative to times)
• times: Specific sample times (alternative to n_samples)

Keyword Arguments

• drive_name: Name for the drive variable (default: :u)
• initial_value: Initial boundary condition (default: pulse(0.0))
• final_value: Final boundary condition (default: pulse(duration))

Example

gaussian = GaussianPulse([1.0, 2.0], 0.1, 1.0)
linear = LinearSplinePulse(gaussian, 50)  # 50 samples

LinearSplinePulse(traj::NamedTrajectory; drive_name=:u)

Construct a LinearSplinePulse from a NamedTrajectory.

Arguments

• traj: NamedTrajectory with control data

Keyword Arguments

• drive_name: Name of the drive component (default: :u)

ZeroOrderPulse{I<:ConstantInterpolation} <: AbstractPulse

Piecewise constant pulse (zero-order hold). The control value at time t is the value at the most recent sample point.

Fields

• controls::I: ConstantInterpolation from DataInterpolations
• duration::Float64: Total pulse duration
• n_drives::Int: Number of control drives
• drive_name::Symbol: Name of the drive variable (default :u)
• initial_value::Vector{Float64}: Initial boundary condition (default: zeros)
• final_value::Vector{Float64}: Final boundary condition (default: zeros)

ZeroOrderPulse(controls::AbstractMatrix, times::AbstractVector; drive_name=:u, initial_value=nothing, final_value=nothing)

Create a zero-order hold pulse from control samples and times.

Arguments

• controls: Matrix of size (n_drives, n_times) with control values
• times: Vector of sample times (must start at 0)

Keyword Arguments

• drive_name: Name of the drive variable (default :u)
• initial_value: Initial boundary condition (default: zeros(n_drives))
• final_value: Final boundary condition (default: zeros(n_drives))

ZeroOrderPulse(traj::NamedTrajectory; drive_name=:u)

Construct a ZeroOrderPulse from a NamedTrajectory.

Arguments

• traj: NamedTrajectory with control data

Keyword Arguments

• drive_name: Name of the drive component (default: :u)

drive_name(pulse::AbstractPulse)

Return the name of the drive variable for this pulse.

duration(pulse::AbstractPulse)

Return the duration of the pulse.

get_knot_count(pulse::AbstractSplinePulse)

Return the number of knots in the spline pulse.

get_knot_derivatives(pulse::CubicSplinePulse)

Return the Hermite tangents at knot points (the du matrix). Only available for CubicSplinePulse.

get_knot_times(pulse::AbstractPulse)

Return the knot/discontinuity times for the pulse.

For piecewise pulses (ZeroOrderPulse, spline pulses), these are the sample times where the control value or its derivative may be discontinuous. For smooth analytic pulses (GaussianPulse, ErfPulse), returns just the endpoints. For CompositePulse, returns the sorted union of all sub-pulse knot times.

get_knot_values(pulse::CubicSplinePulse)

Return the control values at knot points (the u matrix).

n_drives(pulse::AbstractPulse)

Return the number of control drives in the pulse.

sample(pulse::AbstractPulse, times::AbstractVector)

Sample the pulse at the given times. Returns a matrix of size (n_drives, length(times)).

sample(pulse::AbstractPulse; n_samples::Int=100)

Sample the pulse uniformly with n_samples points. Returns (controls, times).

_update_system!(qtraj, sys::QuantumSystem)

Internal method to update the system field in a quantum trajectory. Extended in quantum_trajectories module for specific trajectory types.

extract_globals(traj::NamedTrajectory, names::Vector{Symbol}=Symbol[])

Extract global variables from trajectory as a NamedTuple for easy access. If names is empty, extracts all global variables.

Example

traj = NamedTrajectory(...; global_data=[0.5, 1.0], global_components=(δ=1:1, Ω=2:2))
g = extract_globals(traj)  # (δ = 0.5, Ω = 1.0)

fidelity(ρ::AbstractMatrix{<:Number}, ρ_goal::AbstractMatrix{<:Number})

Calculate the fidelity between two density matrices ρ and ρ_goal.

fidelity(ψ::AbstractVector{<:Number}, ψ_goal::AbstractVector{<:Number})

Calculate the fidelity between two quantum states ψ and ψ_goal.

rollout(qtraj, args...; kwargs...)

Roll out a quantum trajectory with new pulse or ODE parameters. Extended in quantum_trajectories module for specific trajectory types.

rollout!(qtraj, args...; kwargs...)

In-place rollout of quantum trajectory with new pulse or ODE parameters. Extended in quantum_trajectories module for specific trajectory types.

unitary_fidelity(U::AbstractMatrix{<:Number}, U_goal::AbstractMatrix{<:Number})

Calculate the fidelity between unitary operators U and U_goal in the subspace.

update_global_params!(qtraj, traj::NamedTrajectory)

Update the global parameters in the quantum trajectory's system with the optimized values from the NamedTrajectory after optimization. Handles immutable QuantumSystem by reconstructing with updated global_params NamedTuple.

AbstractPiccoloOperator

Union type for operators.

EmbeddedOperator

Embedded operator type to represent an operator embedded in a subspace of a larger quantum system.

Fields

• operator::Matrix{<:Number}: Embedded operator of size prod(subsystem_levels) x prod(subsystem_levels).
• subspace::Vector{Int}: Indices of the subspace the operator is embedded in.
• subsystem_levels::Vector{Int}: Levels of the subsystems in the composite system.

EmbeddedOperator(
    subspace_operator::AbstractMatrix{<:Number},
    subsystem_indices::AbstractVector{Int},
    subspaces::AbstractVector{<:AbstractVector{Int}},
    subsystem_levels::AbstractVector{Int}
)

Embed the subspace_operator into the provided subspaces of a composite system, where the subsystem_indices list the subspaces at which the operator is defined, and the subsystem_levels list the levels of the subsystems in which the operator is embedded.

EmbeddedOperator(
    subspace_operator::AbstractMatrix{<:Number},
    subsystem_indices::AbstractVector{Int},
    subspaces::AbstractVector{<:AbstractVector{Int}},
    composite_system::CompositeQuantumSystem
)

Embed the subspace_operator into the provided subspaces of a composite system.

EmbeddedOperator(subspace_operator::AbstractMatrix{<:Number}, system::QuantumSystem; kwargs...)

Embed the subspace_operator into a quantum system.

EmbeddedOperator(subspace_operator::Matrix{<:Number}, subspace::AbstractVector{Int}, subsystem_levels::AbstractVector{Int})

Create an embedded operator. The operator is embedded at the subspace of the system spanned by the subsystem_levels.

embed(subspace_operator::AbstractMatrix{<:Number}, embedded_operator::EmbeddedOperator)

Embed the subspace_operator in the subspace of a larger embedded_operator.

embed(operator::AbstractMatrix{<:Number}, subspace::AbstractVector{Int}, levels::Int)

Embed an operator in the subspace of a larger matrix of size levels x levels.

get_enr_subspace_indices(excitation_restriction::Int, subsystem_levels::AbstractVector{Int})

Get the indices for the subspace of the quantum system with an excitation restriction.

get_iso_vec_leakage_indices(subspace::AbstractVector{Int}, levels::Int)
get_iso_vec_leakage_indices(subspaces::AbstractVector{<:AbstractVector{Int}}, subsystem_levels::AbstractVector{Int})
get_iso_vec_leakage_indices(subsystem_levels::AbstractVector{Int}; subspace=1:2)
get_iso_vec_leakage_indices(op::EmbeddedOperator)

Get the indices for the leakage in the isomorphic vector space for operators.

get_iso_vec_subspace_indices(subspace::AbstractVector{Int}, subsystem_levels::AbstractVector{Int})
get_iso_vec_subspace_indices(op::EmbeddedOperator)

Get the indices for the subspace in the isomorphic vector space for operators.

get_leakage_indices(subspace::AbstractVector{Int}, levels::Int)
get_leakage_indices(subspaces::AbstractVector{<:AbstractVector{Int}}, subsystem_levels::AbstractVector{Int})
get_leakage_indices(subsystem_levels::AbstractVector{Int}; subspace=1:2)
get_leakage_indices(op::EmbeddedOperator)

Get the indices for the states that are outside of the provided subspace of the quantum system.

get_subspace_indices(subspace::AbstractVector{Int}, levels::Int)
get_subspace_indices(subspaces::Vector{<:AbstractVector{Int}}, subsystem_levels::AbstractVector{Int})
get_subspace_indices(subsystem_levels::AbstractVector{Int}; subspace=1:2)
get_subspace_indices(op::EmbeddedOperator)

Get the indices for the provided subspace of the quantum system.

unembed(matrix::AbstractMatrix{<:Number}, subspace::AbstractVector{Int})

Unembed a subspace operator from the matrix. This is equivalent to calling matrix[subspace, subspace].

unembed(op::AbstractMatrix, embedded_op::EmbeddedOperator)

Unembed a sub-matrix from the op at the subspace defined by embedded_op.

unembed(embedded_op::EmbeddedOperator)

Unembed an embedded operator, returning the original operator.

operator_to_iso_vec(op::EmbeddedOperator)

Convert an EmbeddedOperator into a real vector by converting its embedded operator matrix.

lift_operator(operator::AbstractMatrix{<:Number}, i::Int, subsystem_levels::Vector{Int})
lift_operator(operator::AbstractMatrix{<:Number}, i::Int, n_qubits::Int; kwargs...)
lift_operator(operators::AbstractVector{<:AbstractMatrix{T}}, indices::AbstractVector{Int}, subsystem_levels::Vector{Int})
lift_operator(operators::AbstractVector{<:AbstractMatrix{T}}, indices::AbstractVector{Int}, n_qubits::Int; kwargs...)
lift_operator(operator::AbstractMatrix{T}, indices::AbstractVector{Int}, subsystem_levels::AbstractVector{Int})
lift_operator(operator::AbstractMatrix{T}, indices::AbstractVector{Int}, n_qubits::Int; kwargs...)

Lift an operator acting on the i-th subsystem within subsystem_levels to an operator acting on the entire system spanning subsystem_levels.

direct_sum(A::AbstractMatrix, B::AbstractMatrix)

Returns the direct sum of two matrices.

direct_sum(Ã⃗::AbstractVector, B̃⃗::AbstractVector)

Returns the direct sum of two iso_vec operators.

direct_sum(sys1::QuantumSystem, sys2::QuantumSystem)

Returns the direct sum of two QuantumSystem objects.

Constructs a new system where the Hilbert space is the direct sum of the two input systems: H = H₁ ⊕ H₂ = [H₁ 0 ] [0 H₂]

Both systems must have the same number of drives. The resulting system uses sys1's drive_bounds.

Example

sys1 = QuantumSystem([PAULIS[:X]], [(-1.0, 1.0)])
sys2 = QuantumSystem([PAULIS[:Y]], [(-1.0, 1.0)])
sys_combined = direct_sum(sys1, sys2)

direct_sum(A::SparseMatrixCSC, B::SparseMatrixCSC)

Returns the direct sum of two sparse matrices.

G(H::AbstractMatrix)::Matrix{Float64}

Returns the isomorphism of $-iH$, i.e. $G(H) = \text{iso}(-iH)$.

See also Isomorphisms.iso, Isomorphisms.H.

H(G::AbstractMatrix{<:Real})

Returns the inverse of $G(H) = iso(-iH)$, i.e. returns H.

See also Isomorphisms.iso, Isomorphisms.G.

ad_vec(H::AbstractMatrix{ℂ}; anti::Bool=false) where ℂ <: Number

Returns the vectorized adjoint action of a matrix H:

\[\text{ad_vec}(H) = \mqty(1 & 0 \\ 0 & 1) \otimes H - (-1)^{\text{anti}} \mqty(0 & 1 \\ 1 & 0) \otimes H^*\]

bloch_to_ket(v::AbstractVector{<:Real}; digits=6)

Convert a Bloch vector to a ket (up to global phase).

compact_iso_to_density(x::AbstractVector{<:Real})

Reconstructs a Hermitian density matrix $ρ$ from its compact real isomorphism vector $x$ of length $n^2$.

See also density_to_compact_iso.

density_lift_matrix(n::Int)

Returns the sparse lift matrix $L \in \mathbb{R}^{2n^2 \times n^2}$ that maps the compact density isomorphism to the full $\text{iso\_vec}$ representation.

Given a compact vector $x = \text{density\_to\_compact\_iso}(ρ)$, the full representation is recovered as $L x = \text{density\_to\_iso\_vec}(ρ)$.

The matrix has $n(2n-1)$ nonzero entries.

See also density_projection_matrix, density_to_compact_iso.

density_projection_matrix(n::Int)

Returns the sparse projection matrix $P \in \mathbb{R}^{n^2 \times 2n^2}$ that maps the full $\text{iso\_vec}$ representation to the compact density isomorphism.

Given a full vector $v = \text{density\_to\_iso\_vec}(ρ)$ for Hermitian $ρ$, the compact representation is $P v = \text{density\_to\_compact\_iso}(ρ)$.

Satisfies $PL = I_{n^2}$ where $L$ = density_lift_matrix.

See also density_lift_matrix, density_to_compact_iso.

density_to_compact_iso(ρ::AbstractMatrix{<:Number})

Returns a compact real isomorphism vector of length $n^2$ for a Hermitian density matrix $ρ$ of size $n \times n$. Exploits Hermiticity ($ρ = ρ^\dagger$) to halve the representation size compared to density_to_iso_vec.

The compact vector is ordered as:

1. $\text{Re}(ρ_{jk})$ for $j \leq k$ (upper triangle, column-major): $n(n+1)/2$ entries
2. $\text{Im}(ρ_{jk})$ for $j < k$ (strict upper triangle, column-major): $n(n-1)/2$ entries

See also compact_iso_to_density, density_lift_matrix, density_projection_matrix.

density_to_iso_vec(ρ::AbstractMatrix{<:Number})

Returns the isomorphism ρ⃗̃ = ket_to_iso(vec(ρ)) of a density matrix ρ

iso(H::AbstractMatrix{<:Number})

Returns the isomorphism of $H$:

\[iso(H) = \widetilde{H} = \mqty(1 & 0 \\ 0 & 1) \otimes \Re(H) + \mqty(0 & -1 \\ 1 & 0) \otimes \Im(H)\]

where $\Im(H)$ and $\Re(H)$ are the imaginary and real parts of $H$ and the tilde indicates the standard isomorphism of a complex valued matrix:

\[\widetilde{H} = \mqty(1 & 0 \\ 0 & 1) \otimes \Re(H) + \mqty(0 & -1 \\ 1 & 0) \otimes \Im(H)\]

See also Isomorphisms.G, Isomorphisms.H.

iso_D(L::AbstractMatrix{ℂ}) where ℂ <: Number

Returns the isomorphic representation of the Lindblad dissipator L.

iso_operator_to_iso_vec(Ũ::AbstractMatrix{ℝ}) where ℝ <: Real

Convert a real matrix Ũ representing an isomorphism operator into a real vector.

iso_operator_to_operator(Ũ)

iso_to_ket(ψ̃::AbstractVector{<:Real})

Convert a real isomorphism vector ψ̃ into a ket vector.

iso_vec_to_density(ρ⃗̃::AbstractVector{<:Real})

Returns the density matrix ρ from its isomorphism ρ⃗̃

iso_vec_to_iso_operator(Ũ⃗::AbstractVector{ℝ}) where ℝ <: Real

Convert a real vector Ũ⃗ into a real matrix representing an isomorphism operator.

iso_vec_to_operator(Ũ⃗::AbstractVector{ℝ}) where ℝ <: Real

Convert a real vector Ũ⃗ into a complex matrix representing an operator.

ket_to_bloch(ψ::AbstractVector{<:Number})

Convert a ket to a Bloch vector representation.

ket_to_iso(ψ::AbstractVector{<:Number})

Convert a ket vector ψ into a complex vector with real and imaginary parts.

mat(x::AbstractVector)

Convert a vector x into a square matrix. The length of x must be a perfect square.

operator_to_iso_operator(U)

operator_to_iso_vec(U::AbstractMatrix{ℂ}) where ℂ <: Number

Convert a complex matrix U representing an operator into a real vector.

var_G(G::AbstractMatrix{<:Real}, G_vars::AbstractVector{<:AbstractMatrix{<:Real}})

Returns the variational generator of G with variational derivatives, G_vars.

The variational generator is

\[\text{var}_G(G, [G_a, G_b]) = \mqty( G & 0 & 0 \\ G_a & G & 0 \\ G_b & 0 & G )\]

where G is the isomorphism of a Hamiltonian and G_a and G_b are the variational derivatives of G for parameters a and b, respectively.

QuantumControlProblem{QT<:AbstractQuantumTrajectory}

Wrapper combining quantum trajectory information with trajectory optimization problem.

This type enables:

• Type-stable dispatch on quantum trajectory type (Unitary, Ket, Density)
• Clean separation of quantum information (system, goal) from optimization details
• Composable problem transformations (e.g., SmoothPulseProblem → MinimumTimeProblem)

Fields

• qtraj::QT: Quantum trajectory containing system, goal, and quantum state information
• prob::DirectTrajOptProblem: Direct trajectory optimization problem with objective, dynamics, constraints

Construction

Typically created via problem templates:

qtraj = UnitaryTrajectory(sys, U_goal, N)
qcp = SmoothPulseProblem(qtraj; Q=100.0, R=1e-2)

Accessors

• get_trajectory(qcp): Get the NamedTrajectory
• get_system(qcp): Get the QuantumSystem
• get_goal(qcp): Get the goal state/unitary
• state_name(qcp): Get the state variable name
• drive_name(qcp): Get the control variable name

Solving

solve!(qcp; max_iter=100, verbose=true)

solve!(qcp::QuantumControlProblem; sync::Bool=true, kwargs...)

Solve the quantum control problem by forwarding to the inner DirectTrajOptProblem.

Arguments

• sync::Bool=true: If true, call sync_trajectory! after solving to update qtraj.trajectory with physical control values. Set to false to skip synchronization (e.g., for debugging).

All other keyword arguments are passed to the DirectTrajOpt solver.

get_trajectory(qcp::QuantumControlProblem)

Get the NamedTrajectory from the optimization problem.

sync_trajectory!(qcp::QuantumControlProblem)

Update the quantum trajectory in-place from the optimized control values.

After optimization, this function:

1. Extracts the optimized controls from prob.trajectory (unadapting if needed)
2. Creates a new pulse with those controls via extract_pulse
3. Re-solves the ODE to get the updated quantum evolution
4. Replaces qtraj with the new quantum trajectory

This gives you access to the continuous-time ODE solution with the optimized controls, allowing you to:

• Evaluate the fidelity via fidelity(qcp.qtraj)
• Sample the quantum state at any time via qcp.qtraj(t)
• Get the optimized pulse via get_pulse(qcp.qtraj)

Example

solve!(qcp; max_iter=100)  # Automatically calls sync_trajectory!
fid = fidelity(qcp.qtraj)  # Evaluate fidelity with continuous-time solution
pulse = get_pulse(qcp.qtraj)  # Get the optimized pulse

drive_name(qcp::QuantumControlProblem)

Get the control variable name from the quantum trajectory.

get_goal(qcp::QuantumControlProblem)

Get the goal state/operator from the quantum trajectory.

get_system(qcp::QuantumControlProblem)

Get the QuantumSystem from the quantum trajectory.

state_name(qcp::QuantumControlProblem)

Get the state variable name from the quantum trajectory.

fidelity(qcp::QuantumControlProblem; kwargs...)

Compute the fidelity of the quantum trajectory.

This is a convenience wrapper that forwards to fidelity(qcp.qtraj; kwargs...).

Example

solve!(qcp)
fid = fidelity(qcp)  # Equivalent to fidelity(qcp.qtraj)

BangBangPulseProblem(qtraj::AbstractQuantumTrajectory{<:ZeroOrderPulse}, N::Int; kwargs...)

Construct a QuantumControlProblem that promotes bang-bang (piecewise-constant, few-switch) controls by penalizing $\|du\|_1$ via an exact slack reformulation.

Unlike SmoothPulseProblem (which uses 2 derivative levels with L2 regularization), this stores only 1 derivative (du) and uses slack variables to impose an exact L1 penalty on it, promoting sparsity in du and thus fewer switches.

L1 penalty via slack variables

Introduces slack variables $s \geq 0$ (same dimension as du) and enforces:

\[|du_{k,i}| \leq s_{k,i}\]

Then minimizes the linear cost $R_{du} \sum_k \sum_i s_{k,i} \Delta t_k$. At optimality, $s = |du|$, giving the exact L1 norm.

Arguments

• qtraj::AbstractQuantumTrajectory{<:ZeroOrderPulse}: Quantum trajectory with piecewise constant pulse
• N::Int: Number of timesteps for discretization

Keyword Arguments

• integrator::Union{Nothing, AbstractIntegrator, Vector{<:AbstractIntegrator}}=nothing: Optional custom integrator(s). If not provided, uses BilinearIntegrator.
• global_names::Union{Nothing, Vector{Symbol}}=nothing: Names of global variables to optimize. Requires a custom integrator.
• global_bounds::Union{Nothing, Dict{Symbol, Union{Float64, Tuple{Float64, Float64}}}}=nothing: Bounds for global variables.
• du_bound::Float64=Inf: Bound on discrete first derivative
• Δt_bounds::Union{Nothing, Tuple{Float64, Float64}}=nothing: Timestep bounds
• Q::Float64=100.0: Weight on infidelity/objective
• R::Float64=1e-2: Default regularization weight
• R_u::Union{Float64, Vector{Float64}}=0.0: L2 weight on control amplitude (defaults to 0 — bang-bang needs no amplitude regularization)
• R_du::Union{Float64, Vector{Float64}}=R: L1 weight on du (applied to slacks)
• constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]: Additional constraints
• piccolo_options::PiccoloOptions=PiccoloOptions(): Piccolo solver options

Returns

• QuantumControlProblem: Wrapper containing quantum trajectory and optimization problem

Examples

sys = QuantumSystem(H_drift, H_drives, drive_bounds)
pulse = ZeroOrderPulse(0.1 * randn(n_drives, N), collect(range(0.0, T, length=N)))
qtraj = UnitaryTrajectory(sys, pulse, U_goal)
qcp = BangBangPulseProblem(qtraj, N; Q=100.0, R_du=1e-1)
solve!(qcp; max_iter=200)

See also: SmoothPulseProblem for smooth (L2-regularized) controls.

BangBangPulseProblem(qtraj::MultiKetTrajectory{<:ZeroOrderPulse}, N::Int; kwargs...)

Construct a QuantumControlProblem for bang-bang pulse optimization over an ensemble of ket state transfers with piecewise constant controls.

See the single-trajectory method for full documentation of keyword arguments.

BangBangPulseProblem(qtraj::AbstractQuantumTrajectory, N::Int; kwargs...)

Fallback method that provides helpful error for non-ZeroOrderPulse types.

MinimumTimeProblem(qcp::QuantumControlProblem; kwargs...)

Convert an existing quantum control problem to minimum-time optimization.

IMPORTANT: This function requires an existing QuantumControlProblem (e.g., from SmoothPulseProblem). It cannot be created directly from a quantum trajectory. The workflow is:

1. Create base problem with SmoothPulseProblem (or similar)
2. Solve base problem to get feasible solution
3. Convert to minimum-time with MinimumTimeProblem

This ensures the problem starts from a good initialization and maintains solution quality through the final fidelity constraint.

Type Dispatch

Automatically handles different quantum trajectory types through the type parameter:

• QuantumControlProblem{UnitaryTrajectory} → Uses FinalUnitaryFidelityConstraint
• QuantumControlProblem{KetTrajectory} → Uses FinalKetFidelityConstraint
• QuantumControlProblem{DensityTrajectory} → Not yet implemented

The optimization problem is:

\[\begin{aligned} \underset{\vec{\tilde{q}}, u, \Delta t}{\text{minimize}} & \quad J_{\text{original}}(\vec{\tilde{q}}, u) + D \sum_t \Delta t_t \\ \text{ subject to } & \quad \text{original dynamics \& constraints} \\ & F_{\text{final}} \geq F_{\text{threshold}} \\ & \quad \Delta t_{\text{min}} \leq \Delta t_t \leq \Delta t_{\text{max}} \\ \end{aligned}\]

where q represents the quantum state (unitary, ket, or density matrix).

Arguments

• qcp::QuantumControlProblem: Existing quantum control problem to convert

Keyword Arguments

• final_fidelity::Float64=0.99: Minimum fidelity constraint at final time
• D::Float64=100.0: Weight on minimum-time objective ∑Δt
• piccolo_options::PiccoloOptions=PiccoloOptions(): Piccolo solver options

Returns

• QuantumControlProblem: New problem with minimum-time objective and fidelity constraint

Examples

# Standard workflow
sys = QuantumSystem(H_drift, H_drives, drive_bounds)
pulse = ZeroOrderPulse(0.1 * randn(n_drives, N), collect(range(0.0, T, length=N)))
qtraj = UnitaryTrajectory(sys, pulse, U_goal)

# Step 1: Create and solve base smooth pulse problem (with Δt_bounds for free time)
qcp_smooth = SmoothPulseProblem(qtraj, N; Q=100.0, R=1e-2, Δt_bounds=(0.01, 0.5))
solve!(qcp_smooth; max_iter=100)

# Step 2: Convert to minimum-time
qcp_mintime = MinimumTimeProblem(qcp_smooth; final_fidelity=0.99, D=100.0)
solve!(qcp_mintime; max_iter=100)

# Compare durations
duration_before = sum(get_timesteps(get_trajectory(qcp_smooth)))
duration_after = sum(get_timesteps(get_trajectory(qcp_mintime)))
@assert duration_after <= duration_before

# Nested transformations also work
qcp_final = MinimumTimeProblem(
    RobustnessProblem(qcp_smooth);  # Future feature
    final_fidelity=0.95
)

Convenience Constructors

You can also update the goal when creating minimum-time problem:

# Different goal for minimum-time optimization
qcp_mintime = MinimumTimeProblem(qcp_smooth; goal=U_goal_new, final_fidelity=0.98)

SamplingProblem(qcp::QuantumControlProblem, systems::Vector{<:AbstractQuantumSystem}; kwargs...)

Construct a SamplingProblem from an existing QuantumControlProblem and a list of systems.

This creates a robust optimization problem where the controls are shared across all systems, but each system evolves according to its own dynamics. The objective is the weighted sum of fidelity objectives for each system.

Arguments

• qcp::QuantumControlProblem: The base problem (defines nominal trajectory, objective, etc.)
• systems::Vector{<:AbstractQuantumSystem}: List of systems to optimize over

Keyword Arguments

• weights::Vector{Float64}=fill(1.0, length(systems)): Weights for each system
• Q::Float64=100.0: Weight on infidelity objective (explicit, not extracted from base problem)
• piccolo_options::PiccoloOptions=PiccoloOptions(): Options for the solver

Returns

• QuantumControlProblem{SamplingTrajectory}: A new problem with the sampling trajectory

SmoothPulseProblem(qtraj::AbstractQuantumTrajectory{<:ZeroOrderPulse}, N::Int; kwargs...)

Construct a QuantumControlProblem for smooth pulse optimization with piecewise constant controls.

Note: This problem template is for ZeroOrderPulse only. For spline-based pulses (LinearSplinePulse, CubicSplinePulse), use SplinePulseProblem instead.

The problem adds discrete derivative variables (du, ddu) that:

• Regularize control changes between timesteps
• Enforce smoothness via DerivativeIntegrator constraints

Arguments

• qtraj::AbstractQuantumTrajectory{<:ZeroOrderPulse}: Quantum trajectory with piecewise constant pulse
• N::Int: Number of timesteps for discretization

Keyword Arguments

• integrator::Union{Nothing, AbstractIntegrator, Vector{<:AbstractIntegrator}}=nothing: Optional custom integrator(s). If not provided, uses BilinearIntegrator (which does not support global variables). A custom integrator is required when global_names is specified.
• global_names::Union{Nothing, Vector{Symbol}}=nothing: Names of global variables to optimize. Requires a custom integrator (e.g., HermitianExponentialIntegrator from Piccolissimo) that supports global variables.
• global_bounds::Union{Nothing, Dict{Symbol, Union{Float64, Tuple{Float64, Float64}}}}=nothing: Bounds for global variables. Keys are variable names, values are either a scalar (symmetric bounds ±value) or a tuple (lower, upper).
• du_bound::Float64=Inf: Bound on discrete first derivative (controls jump rate)
• ddu_bound::Float64=1.0: Bound on discrete second derivative (controls acceleration)
• Q::Float64=100.0: Weight on infidelity/objective
• R::Float64=1e-2: Weight on regularization terms (u, u̇, ü)
• R_u::Union{Float64, Vector{Float64}}=R: Weight on control regularization
• R_du::Union{Float64, Vector{Float64}}=R: Weight on first derivative regularization
• R_ddu::Union{Float64, Vector{Float64}}=R: Weight on second derivative regularization
• constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]: Additional constraints
• piccolo_options::PiccoloOptions=PiccoloOptions(): Piccolo solver options

Returns

• QuantumControlProblem: Wrapper containing quantum trajectory and optimization problem

Examples

# Unitary gate synthesis with piecewise constant pulse
sys = QuantumSystem(H_drift, H_drives, drive_bounds)
pulse = ZeroOrderPulse(0.1 * randn(n_drives, N), collect(range(0.0, T, length=N)))
qtraj = UnitaryTrajectory(sys, pulse, U_goal)
qcp = SmoothPulseProblem(qtraj, N; Q=100.0, R=1e-2)
solve!(qcp; max_iter=100)

# Quantum state transfer
pulse = ZeroOrderPulse(0.1 * randn(n_drives, N), collect(range(0.0, T, length=N)))
qtraj = KetTrajectory(sys, pulse, ψ_init, ψ_goal)
qcp = SmoothPulseProblem(qtraj, N; Q=50.0, R=1e-3)
solve!(qcp)

See also: SplinePulseProblem for spline-based pulses.

SmoothPulseProblem(qtraj::MultiKetTrajectory{<:ZeroOrderPulse}, N::Int; kwargs...)

Construct a QuantumControlProblem for smooth pulse optimization over an ensemble of ket state transfers with piecewise constant controls.

This handles the case where you want to optimize a single pulse that achieves multiple state transfers simultaneously (e.g., |0⟩→|1⟩ and |1⟩→|0⟩ for an X gate via state transfer).

Note: This problem template is for ZeroOrderPulse only. For spline-based pulses, use SplinePulseProblem instead.

Arguments

• qtraj::MultiKetTrajectory{<:ZeroOrderPulse}: Ensemble of ket state transfers with piecewise constant pulse
• N::Int: Number of timesteps for the discretization

Keyword Arguments

• integrator::Union{Nothing, AbstractIntegrator, Vector{<:AbstractIntegrator}}=nothing: Optional custom integrator(s). If not provided, the default BilinearIntegrator is used. When global_names is specified, you must supply a custom integrator here (i.e., do not rely on the default BilinearIntegrator) that supports global variables.
• global_names::Union{Nothing, Vector{Symbol}}=nothing: Names of global variables to optimize. Requires a custom integrator provided via integrator (e.g., HermitianExponentialIntegrator from Piccolissimo) that supports global variables.
• global_bounds::Union{Nothing, Dict{Symbol, Union{Float64, Tuple{Float64, Float64}}}}=nothing: Bounds for global variables. Keys are variable names, values are either a scalar (symmetric bounds ±value) or a tuple (lower, upper).
• du_bound::Float64=Inf: Bound on discrete first derivative
• ddu_bound::Float64=1.0: Bound on discrete second derivative
• Q::Float64=100.0: Weight on infidelity/objective
• R::Float64=1e-2: Weight on regularization terms (u, u̇, ü)
• R_u::Union{Float64, Vector{Float64}}=R: Weight on control regularization
• R_du::Union{Float64, Vector{Float64}}=R: Weight on first derivative regularization
• R_ddu::Union{Float64, Vector{Float64}}=R: Weight on second derivative regularization
• constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]: Additional constraints
• piccolo_options::PiccoloOptions=PiccoloOptions(): Piccolo solver options

Returns

• QuantumControlProblem{MultiKetTrajectory}: Wrapper containing ensemble trajectory and optimization problem

Examples

# Create ensemble for X gate via state transfer
sys = QuantumSystem(H_drift, H_drives, drive_bounds)
pulse = ZeroOrderPulse(0.1 * randn(n_drives, N), collect(range(0.0, T, length=N)))

ψ0 = ComplexF64[1.0, 0.0]
ψ1 = ComplexF64[0.0, 1.0]

ensemble_qtraj = MultiKetTrajectory(sys, pulse, [ψ0, ψ1], [ψ1, ψ0])
qcp = SmoothPulseProblem(ensemble_qtraj, N; Q=100.0, R=1e-2)
solve!(qcp; max_iter=100)

See also: SplinePulseProblem for spline-based pulses.

SmoothPulseProblem(qtraj::AbstractQuantumTrajectory, N::Int; kwargs...)

Fallback method that provides helpful error for non-ZeroOrderPulse types.

SplinePulseProblem(qtraj::MultiKetTrajectory{<:AbstractSplinePulse}; kwargs...)
SplinePulseProblem(qtraj::MultiKetTrajectory{<:AbstractSplinePulse}, N::Int; kwargs...)
SplinePulseProblem(qtraj::MultiKetTrajectory{<:AbstractSplinePulse}, times::AbstractVector; kwargs...)

Create a spline-based trajectory optimization problem for ensemble ket state transfers.

Uses coherent fidelity objective (phases must align) for gate implementation.

Arguments

• qtraj::MultiKetTrajectory{<:AbstractSplinePulse}: Ensemble trajectory with spline pulse
• N_or_times: One of:
  • nothing (default): Use native knot times from spline pulse
  • N::Int: Number of uniformly spaced timesteps
  • times::AbstractVector: Specific sample times

Keyword Arguments

Same as the base SplinePulseProblem method.

SplinePulseProblem(qtraj::AbstractQuantumTrajectory{<:AbstractSplinePulse}; kwargs...)
SplinePulseProblem(qtraj::AbstractQuantumTrajectory{<:AbstractSplinePulse}, N::Int; kwargs...)
SplinePulseProblem(qtraj::AbstractQuantumTrajectory{<:AbstractSplinePulse}, times::AbstractVector; kwargs...)

Construct a QuantumControlProblem for spline-based pulse optimization.

Unlike SmoothPulseProblem (which uses piecewise constant controls with discrete smoothing variables), this problem template is designed for spline-based pulses where the derivative variables (du) are the actual spline coefficients or slopes.

Pulse Type Semantics

LinearSplinePulse: The du variable represents the slope between knots. A DerivativeIntegrator constraint enforces du[k] = (u[k+1] - u[k]) / Δt, making the slopes consistent with the linear interpolation. This constraint ensures mathematical rigor while allowing slope regularization/bounds.

CubicSplinePulse (Hermite spline): The du variable is the tangent/derivative at each knot point, which is a true independent degree of freedom in Hermite interpolation. No DerivativeIntegrator is added - the optimizer can adjust both :u and :du independently.

Mathematical Notes

• LinearSplinePulse: Always adds :du and DerivativeIntegrator to enforce slope consistency
• CubicSplinePulse: :du values are Hermite tangents (unconstrained, only regularized)

Both pulse types always have :du components in the trajectory, simplifying integrator implementations.

Arguments

• qtraj::AbstractQuantumTrajectory{<:AbstractSplinePulse}: Quantum trajectory with spline pulse
• N_or_times: One of:
  • nothing (default): Use native knot times from spline pulse (ideal for warm-starting)
  • N::Int: Number of uniformly spaced timesteps
  • times::AbstractVector: Specific sample times

Keyword Arguments

• integrator::Union{Nothing, AbstractIntegrator, Vector{<:AbstractIntegrator}}=nothing: Optional custom integrator(s). If not provided, uses BilinearIntegrator (which does not support global variables). A custom integrator is required when global_names is specified.
• global_names::Union{Nothing, Vector{Symbol}}=nothing: Names of global variables to optimize. Requires a custom integrator (e.g., SplineIntegrator from Piccolissimo) that supports global variables.
• global_bounds::Union{Nothing, Dict{Symbol, Union{Float64, Tuple{Float64, Float64}}}}=nothing: Bounds for global variables. Keys are variable names, values are either a scalar (symmetric bounds ±value) or a tuple (lower, upper).
• du_bound::Float64=Inf: Bound on derivative (slope) magnitude
• Q::Float64=100.0: Weight on infidelity/objective
• R::Float64=1e-2: Weight on regularization terms
• R_u::Union{Float64, Vector{Float64}}=R: Weight on control regularization
• R_du::Union{Float64, Vector{Float64}}=R: Weight on derivative regularization
• constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]: Additional constraints
• piccolo_options::PiccoloOptions=PiccoloOptions(): Piccolo solver options

Returns

• QuantumControlProblem{<:AbstractQuantumTrajectory}: Wrapper containing trajectory and optimization problem

Examples

# Create system and initial pulse
sys = QuantumSystem(H_drift, H_drives, drive_bounds)
pulse = CubicSplinePulse(u_init, du_init, times)
qtraj = UnitaryTrajectory(sys, pulse, U_goal)

# Use native knot structure (best for warm-starting from saved pulse)
qcp = SplinePulseProblem(qtraj; Q=100.0, du_bound=10.0)

# Or resample to different number of knots
qcp = SplinePulseProblem(qtraj, 50; Q=100.0, du_bound=10.0)

solve!(qcp; max_iter=100)

See also: SmoothPulseProblem for piecewise constant pulses with discrete smoothing.

SplinePulseProblem(qtraj::AbstractQuantumTrajectory, N_or_times; kwargs...)

Fallback method that provides helpful error for non-spline pulse types.

_ensemble_ket_objective(qtraj::MultiKetTrajectory, traj, state_names, weights, goals, Q)

Create a coherent fidelity objective for ensemble state transfers.

For ensemble trajectories (implementing a gate via multiple state transfers), we use coherent fidelity: Fcoherent = |1/n ∑ᵢ ⟨ψᵢgoal|ψᵢ⟩|²

This requires all state overlaps to have aligned phases, which is essential for gate implementation (the gate should have a single global phase).

_final_fidelity_constraint(qtraj::MultiKetTrajectory, final_fidelity, traj)

Create a coherent fidelity constraint for an MultiKetTrajectory.

Uses coherent fidelity: F = |1/n ∑ᵢ ⟨ψᵢ_goal|ψᵢ⟩|²

This enforces that all state transfers have aligned global phases, which is essential when implementing a gate via state transfer (e.g., X gate via |0⟩→|1⟩ and |1⟩→|0⟩).

_make_free_phase_goal(op::EmbeddedOperator)

Build a function θ -> EmbeddedOperator that applies single-qubit Z-phase rotations to the goal gate. For an N-qubit gate, θ has N elements (one phase per qubit).

The phase-adjusted gate is (Z(θ₁) ⊗ Z(θ₂) ⊗ ⋯) ⋅ U_goal, where Z(θ) = diag(1, e^{iθ}).

_make_free_phase_ket_goal(goal, subsystem_levels)

Single-ket version of _make_free_phase_ket_goals. Returns θ -> Vector{ComplexF64}.

_make_free_phase_ket_goals(goals, subsystem_levels)

Build a function θ -> Vector{Vector{ComplexF64}} that applies single-qubit Z-phase rotations to ket goal states in the full Hilbert space.

For an N-qubit system with levels [l₁, l₂, ...], each basis element |s₁,s₂,...,sₙ⟩ gets phase exp(i·Σⱼ θⱼ·δ(sⱼ,1)) — only the |1⟩ level of each qubit acquires a phase.

add_global_bounds_constraints!(constraints, global_bounds, traj; verbose=false)

Add GlobalBoundsConstraint entries for each global variable specified in global_bounds.

Converts bounds from user-friendly formats to the format expected by GlobalBoundsConstraint:

• Float64: Symmetric scalar bounds (applied symmetrically to all dimensions)
• Tuple{Float64, Float64}: Asymmetric scalar bounds (expanded to vectors)
• Vector or Tuple{Vector, Vector}: Already in correct format (passed through)

Modifies constraints in place.

extract_regularization(objective, state_sym::Symbol, new_traj::NamedTrajectory) -> AbstractObjective

Extract regularization terms (non-state-dependent objectives) from a composite objective, filtering to only include terms for variables that exist in the new trajectory.

Used by SamplingProblem to extract shared regularizers (e.g., control penalty) from the base problem while excluding regularizers for variables that don't exist in the sampling trajectory (e.g., :du, :ddu which are added by SmoothPulseProblem).

sampling_state_objective(qtraj, traj, state_sym, Q)

Create the state-dependent objective for a sampling member. Dispatches on quantum trajectory type.

setup_free_phase_globals!(n_qubits, global_data, global_bounds; verbose=false)

Add per-qubit Z-phase variables (φ_1, φ_2, …) to global_data and global_bounds dicts. Returns the vector of phase variable names and the (possibly initialized) dicts.

Mutates global_data and global_bounds in place if they are not nothing; otherwise creates new dicts.

CoherentKetFreePhaseInfidelityObjective(goals_fn, ψ̃_names, θ_names, traj; Q=100.0)

Coherent ket infidelity with optimizable single-qubit Z-phase rotations.

goals_fn(θ) returns phase-rotated goal kets. Phase variables θ are stored in the trajectory's global_data and optimized alongside the pulse.

The objective minimizes: 1 - |1/n Σᵢ ⟨goali(θ)|ψi⟩|²

where goal_i(θ) = Φ(θ) · goal_i with Φ(θ) = Z₁(θ₁) ⊗ Z₂(θ₂) ⊗ ⋯.

Arguments

• goals_fn::Function: Maps phase vector θ to phased goal kets
• ψ̃_names::Vector{Symbol}: Names of isomorphic state variables in trajectory
• θ_names::AbstractVector{Symbol}: Names of phase global variables
• traj::NamedTrajectory: The trajectory

Keyword Arguments

• Q::Float64=100.0: Weight on the infidelity objective

CoherentKetInfidelityObjective(ψ_goals, ψ̃_names, traj; Q=100.0)

Create a terminal objective for coherent ket state infidelity across multiple states.

Coherent fidelity is defined as: Fcoherent = |1/n ∑ᵢ ⟨ψᵢgoal|ψᵢ⟩|²

Unlike incoherent fidelity (average of individual |⟨ψᵢ_goal|ψᵢ⟩|²), coherent fidelity requires all state overlaps to have aligned phases. This is essential when implementing a gate via multiple state transfers - the gate should have a single global phase, not independent phases per state.

Arguments

• ψ_goals::Vector{<:AbstractVector{<:Complex}}: Target ket states
• ψ̃_names::Vector{Symbol}: Names of isomorphic state variables in trajectory
• traj::NamedTrajectory: The trajectory

Keyword Arguments

• Q::Float64=100.0: Weight on the infidelity objective

Example

# For implementing X gate via |0⟩→|1⟩ and |1⟩→|0⟩
goals = [ComplexF64[0, 1], ComplexF64[1, 0]]
names = [:ψ̃1, :ψ̃2]
obj = CoherentKetInfidelityObjective(goals, names, traj; Q=100.0)

DensityMatrixInfidelityObjective(ρ̃_name, ρ_goal, traj; Q=100.0)

Terminal objective for density matrix fidelity using the compact isomorphism.

Minimizes |1 - tr(ρ * ρ_goal)| where ρ is reconstructed from the compact iso vector via compact_iso_to_density.

KetFreePhaseInfidelityObjective(goal_fn, ψ̃_name, θ_names, traj; Q=100.0)

Single-ket infidelity with optimizable Z-phase rotations.

goal_fn(θ) returns the phase-rotated goal ket. Minimizes 1 - |⟨goal(θ)|ψ⟩|².

Arguments

• goal_fn::Function: Maps phase vector θ to phased goal ket
• ψ̃_name::Symbol: Name of isomorphic state variable in trajectory
• θ_names::AbstractVector{Symbol}: Names of phase global variables
• traj::NamedTrajectory: The trajectory

Keyword Arguments

• Q::Float64=100.0: Weight on the infidelity objective

KetInfidelityObjective(ψ_goal, ψ̃_name, traj; Q=100.0)

Create a terminal objective for ket state infidelity with an explicit goal state.

This variant is useful for SamplingProblem and EnsembleTrajectory where the goal is shared across multiple state variables that don't have individual goals in traj.goal.

Arguments

• ψ_goal::AbstractVector{<:Complex}: The target ket state (complex vector)
• ψ̃_name::Symbol: Name of the isomorphic state variable in the trajectory
• traj::NamedTrajectory: The trajectory

Keyword Arguments

• Q::Float64=100.0: Weight on the infidelity objective

KetInfidelityObjective(ψ̃_name, traj; Q=100.0)

Create a terminal objective for ket state infidelity, using the goal from traj.goal[ψ̃_name].

LeakageObjective(indices, name, traj::NamedTrajectory)

Construct a KnotPointObjective that penalizes leakage of name at the knot points specified by times at any indices that are outside the computational subspace.

coherent_ket_fidelity(ψ̃s, ψ_goals)

Compute coherent fidelity across multiple ket states:

F_coherent = |1/n ∑ᵢ ⟨ψᵢ_goal|ψᵢ⟩|²

This requires all overlaps to have consistent phases (global phase alignment), which is necessary for implementing gates via state transfer.

Arguments

• ψ̃s::Vector{<:AbstractVector}: List of isomorphic state vectors
• ψ_goals::Vector{<:AbstractVector{<:Complex}}: List of goal states

FinalCoherentKetFidelityConstraint(ψ_goals, ψ̃_names, final_fidelity, traj)

Create a final fidelity constraint using coherent ket fidelity across multiple states.

Coherent fidelity: F = |1/n ∑ᵢ ⟨ψᵢ_goal|ψᵢ⟩|²

This constraint enforces that all state overlaps have aligned phases, which is essential when implementing a gate via multiple state transfers (e.g., MultiKetTrajectory).

Arguments

• ψ_goals::Vector{<:AbstractVector{<:Complex}}: Target ket states
• ψ̃_names::Vector{Symbol}: Names of isomorphic state variables in trajectory
• final_fidelity::Float64: Minimum fidelity threshold (constraint: F ≥ final_fidelity)
• traj::NamedTrajectory: The trajectory

Example

# For implementing X gate via |0⟩→|1⟩ and |1⟩→|0⟩
goals = [ComplexF64[0, 1], ComplexF64[1, 0]]
names = [:ψ̃1, :ψ̃2]
constraint = FinalCoherentKetFidelityConstraint(goals, names, 0.99, traj)

LeakageConstraint(value, indices, name, traj::NamedTrajectory)

Construct a KnotPointConstraint that bounds leakage of name at the knot points specified by times at any indices that are outside the computational subspace.

PiccoloOptions

Options for the Piccolo quantum optimal control library.

Fields

• verbose::Bool = true: Print verbose output
• timesteps_all_equal::Bool = true: Use equal timesteps
• rollout_integrator::Function = expv: Integrator to use for rollout
• geodesic = true: Use the geodesic to initialize the optimization.
• zero_initial_and_final_derivative::Bool=false: Zero the initial and final control pulse derivatives.
• complex_control_norm_constraint_name::Union{Nothing, Symbol} = nothing: Name of the complex control norm constraint.
• complex_control_norm_constraint_radius::Float64 = 1.0: Radius of the complex control norm constraint.
• bound_state::Bool = false: Bound the state variables <= 1.0.
• leakage_constraint::Bool = false: Suppress leakage with constraint and cost.
• leakage_constraint_value::Float64 = 1e-2: Value for the leakage constraint.
• leakage_cost::Float64 = 1e-2: Leakage suppression parameter.

Calculate Weyl chamber coordinates

plot_state_populations(
    traj::NamedTrajectory;
    state_name::Symbol=:ψ̃,
    state_indices::Union{Nothing, AbstractVector{Int}}=nothing,
    control_name::Symbol=:u,
    subspace::Union{Nothing, AbstractVector{Int}}=nothing,
    kwargs...
)

Plot populations for multiple quantum states stored in a trajectory.

This function visualizes the time evolution of quantum state populations for trajectories containing multiple state trajectories (e.g., from sampling problems where multiple initial states are evolved). States are identified by a common prefix and numeric suffix pattern (e.g., :ψ̃1_system_1, :ψ̃2_system_1, etc.).

Mathematical Background

For a quantum state $|\psi(t)\rangle \in \mathcal{H}$ evolving under the Schrödinger equation, the population in computational basis state $|i\rangle$ is given by

\[P_i(t) = |\langle i|\psi(t)\rangle|^2 = |\psi_i(t)|^2\]

where $\psi_i(t)$ is the $i$-th component of the state vector in the computational basis.

For a normalized state, we have the conservation property:

\[\sum_{i=1}^{n} P_i(t) = \langle\psi(t)|\psi(t)\rangle = 1\]

When multiple states are being optimized simultaneously (as in sampling problems), this function plots the populations for each state, allowing comparison of how different initial conditions evolve under the same control protocol.

Key Properties:

• Populations are real and bounded: $P_i(t) \in [0,1]$
• Total probability is conserved: $\sum_i P_i(t) = 1$
• Can optionally restrict to a subspace (e.g., computational subspace excluding leakage states)

Arguments

• traj::NamedTrajectory: A trajectory containing one or more quantum states in isomorphism representation.

Keyword Arguments

• state_name::Symbol: The base name for state components. The function will find all trajectory components matching this prefix (e.g., :ψ̃ matches :ψ̃1_system_1, :ψ̃2_system_1, etc.). Default is :ψ̃.
• state_indices::Union{Nothing, AbstractVector{Int}}: If provided, only plot states with these indices (e.g., [1, 3] plots only the 1st and 3rd states). If nothing, plots all states matching the prefix. Default is nothing.
• control_name::Symbol: The name of the control signal component to include in the plot. Default is :u.
• subspace::Union{Nothing, AbstractVector{Int}}: If provided, only plot populations for these basis states (e.g., 1:2 for a qubit subspace). Useful for excluding leakage levels. Default is nothing (plot all levels).
• kwargs...: Additional keyword arguments passed to NamedTrajectories.plot.

Returns

• A Makie Figure object containing the population plots for all selected states.

Example

using NamedTrajectories
using Piccolo

# Example: Two initial states evolving under the same controls
N = 100
Δt = 0.1

# Initial states
ψ1_init = ComplexF64[1, 0, 0]  # |0⟩
ψ2_init = ComplexF64[0, 1, 0]  # |1⟩

# Create trajectory with multiple states
traj = NamedTrajectory(
    (
        ψ̃1_system_1 = hcat([ket_to_iso(ψ1_init) for _ in 1:N]...),
        ψ̃2_system_1 = hcat([ket_to_iso(ψ2_init) for _ in 1:N]...),
        u = randn(2, N),
        Δt = fill(Δt, N),
    );
    controls = :u,
    timestep = :Δt,
)

# Plot populations for all states
plot_state_populations(traj)

# Plot only computational subspace (excluding 3rd level)
plot_state_populations(traj; subspace=1:2)

# Plot only first state
plot_state_populations(traj; state_indices=[1])

See also: plot_unitary_populations, NamedTrajectories.plot

plot_unitary_populations(
    traj::NamedTrajectory;
    unitary_columns::AbstractVector{Int}=1:2,
    unitary_name::Symbol=:Ũ⃗,
    control_name::Symbol=:u,
    kwargs...
)

Plot the state populations for specified columns of a unitary operator trajectory.

This function visualizes how the populations (squared magnitudes) of quantum states evolve over time for selected columns of a unitary matrix stored in a NamedTrajectory.

Mathematical Background

For a unitary operator $U(t) \in \mathcal{U}(n)$ evolving under a time-dependent Hamiltonian, this function plots the populations

\[P_{i,j}(t) = |U_{i,j}(t)|^2\]

where $U_{i,j}(t)$ is the $(i,j)$-th element of the unitary matrix at time $t$.

For a quantum system evolving according to the Schrödinger equation

\[\frac{d}{dt}U(t) = -iH(t)U(t), \quad U(0) = I\]

each column $j$ of $U(t)$ represents the time evolution of the initial basis state $|j\rangle$:

\[|\psi_j(t)\rangle = U(t)|j\rangle = \sum_{i=1}^{n} U_{i,j}(t)|i\rangle\]

The population $P_{i,j}(t) = |U_{i,j}(t)|^2$ gives the probability of finding the system in state $|i\rangle$ at time $t$ given that it started in state $|j\rangle$.

Key Properties:

• Unitarity ensures $\sum_{i=1}^{n} P_{i,j}(t) = 1$ for all $j$ and $t$ (probability conservation)
• $P_{i,j}(0) = \delta_{ij}$ (initially in definite state)
• Populations are real and bounded: $P_{i,j}(t) \in [0,1]$

The trajectory stores $U(t)$ in isomorphism representation $\tilde{U}(t)$, a vectorized form that preserves the operator structure while enabling efficient optimization algorithms.

Arguments

• traj::NamedTrajectory: A trajectory containing a unitary operator in isomorphism representation.

Keyword Arguments

• unitary_columns::AbstractVector{Int}: Indices of unitary matrix columns to plot. Each column $j$ corresponds to the evolution of basis state $|j\rangle$. Default is 1:2.
• unitary_name::Symbol: The name of the unitary operator component in the trajectory, stored as an isomorphism vector ($\tilde{U}$). Default is :Ũ⃗.
• control_name::Symbol: The name of the control signal component to include in the plot, typically the time-dependent control parameters $u(t)$ in $H(t) = H_0 + \sum_k u_k(t) H_k$. Default is :u.
• kwargs...: Additional keyword arguments passed to NamedTrajectories.plot, such as xlims, ylims, or Makie-specific plotting options.

Returns

• A Makie Figure object containing the population plots.

Example

using NamedTrajectories
using Piccolo
using Piccolo

# Define Hamiltonian: H = X + a₁(t)Z + a₂(t)Y
H_drift = PAULIS[:X]
H_drives = [PAULIS[:Z], PAULIS[:Y]]

# Generate control trajectory
N, T = 100, 10.0
times = collect(range(0, T, length = N))
u = 0.1 * randn(2, N)

# Generate unitaries
Us = exp.(-im * [(H_drift + sum(u[:, k] .* H_drives)) * times[k] for k = 1:N])

# Create trajectory
traj = NamedTrajectory(
    (
        Ũ⃗ = hcat(operator_to_iso_vec.(Us)...),
        u = u,
        Δt = fill(T / N, N),
    );
    controls = :u,
    timestep = :Δt,
)

# Plot populations for first two columns
plot_unitary_populations(traj)

# Plot only the first column
plot_unitary_populations(traj; unitary_columns=[1])

See also: NamedTrajectories.plot

plot the weyl chamber and record animated trajectory