Library Reference <a class="docs-heading-anchor-permalink" href="#Library-Reference" title="Permalink">

Piccolo.Quantum.QuantumSystems.CompositeQuantumSystem — Type

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

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[])

Piccolo.Quantum.QuantumSystems.CompositeQuantumSystem — Method

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[])

Piccolo.Quantum.QuantumSystems.CompositeQuantumSystem — Method

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[])

Piccolo.Quantum.QuantumSystems.CompositeQuantumSystem — Method

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[])

Piccolo.Quantum.QuantumSystems.CompositeQuantumSystem — Method

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

Piccolo.Quantum.QuantumSystems.DriveBounds — Type

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)]

Piccolo.Quantum.QuantumSystems.OpenQuantumSystem — Type

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{SparseMatrixCSC{ComplexF64, Int}}: The drive Hamiltonians
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

Quantum System Templates

Piccolo.Quantum.QuantumSystemTemplates.CatSystem — Method

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

Returns an OpenQuantumSystem for a quantum cat.

Piccolo.Quantum.QuantumSystemTemplates.IonChainSystem — Method

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).

Piccolo.Quantum.QuantumSystemTemplates.MolmerSorensenCoupling — Method

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.

Piccolo.Quantum.QuantumSystemTemplates.MultiTransmonSystem — Method

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.

Piccolo.Quantum.QuantumSystemTemplates.RadialMSGateSystem — Method

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

Motional closure: All 2N modes must satisfy $|\alpha_k(\tau)| \approx 0$
Target mode: Choose one mode (e.g., k=1) as primary entangling mode
Spectator modes: Other modes should remain minimally excited
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)

Piccolo.Quantum.QuantumSystemTemplates.RadialMSGateSystemWithPhase — Method

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)

Piccolo.Quantum.QuantumSystemTemplates.RydbergChainSystem — Method

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.

Piccolo.Quantum.QuantumSystemTemplates.TransmonCavitySystem — Method

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:

$\hat{a}^\dagger + \hat{a}$ - Real transmon drive
$i(\hat{a}^\dagger - \hat{a})$ - Imaginary transmon drive
$\hat{b}^\dagger + \hat{b}$ - Real cavity drive
$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)

Piccolo.Quantum.QuantumSystemTemplates.TransmonDipoleCoupling — Function

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.

Piccolo.Quantum.QuantumSystemTemplates.TransmonSystem — Method

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.

Piccolo.Quantum.QuantumSystemTemplates.lift — Method

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).

Quantum System Utilities

Piccolo.Quantum.QuantumSystemUtils.is_linearly_dependent — Method

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

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

Piccolo.Quantum.QuantumSystemUtils.is_reachable — Method

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

Piccolo.Quantum.QuantumSystemUtils.is_reachable — Method

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

Piccolo.Quantum.QuantumSystemUtils.operator_algebra — Method

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

Gates and Pauli Matrices

Piccolo.Quantum.Gates.GATES — Constant

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.

Piccolo.Quantum.Gates.PAULIS — Constant

The 2×2 Pauli matrics and identity.

Quantum Object Utilities

Piccolo.Quantum.QuantumObjectUtils.annihilate — Method

annihilate(levels::Int)

Get the annihilation operator for a system with levels.

Piccolo.Quantum.QuantumObjectUtils.create — Method

create(levels::Int)

Get the creation operator for a system with levels.

Piccolo.Quantum.QuantumObjectUtils.haar_identity — Method

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.

Piccolo.Quantum.QuantumObjectUtils.haar_random — Method

haar_random(n::Int)

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

Piccolo.Quantum.QuantumObjectUtils.ket_from_bitstring — Method

ket_from_bitstring(ket::String)

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

Piccolo.Quantum.QuantumObjectUtils.ket_from_string — Method

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.

Piccolo.Quantum.QuantumObjectUtils.operator_from_string — Method

operator_from_string(operator::String; lookup=PAULIS)

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

Trajectories

NamedTrajectories.StructNamedTrajectory.NamedTrajectory — Type

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.

NamedTrajectories.StructNamedTrajectory.NamedTrajectory — Type

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

Convert a DensityTrajectory to a NamedTrajectory for optimization.

Stored Variables

ρ⃗̃: Vectorized isomorphism of the density matrix
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).

NamedTrajectories.StructNamedTrajectory.NamedTrajectory — Type

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.

NamedTrajectories.StructNamedTrajectory.NamedTrajectory — Type

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.

NamedTrajectories.StructNamedTrajectory.NamedTrajectory — Method

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).

Piccolo.Quantum.QuantumTrajectories.AbstractQuantumTrajectory — Type

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)

Piccolo.Quantum.QuantumTrajectories.DensityTrajectory — Type

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.

Piccolo.Quantum.QuantumTrajectories.DensityTrajectory — Method

DensityTrajectory(system, initial, goal, T::Real; drive_name=:u, algorithm=Tsit5())

Convenience constructor that creates a zero pulse of duration T.

Piccolo.Quantum.QuantumTrajectories.DensityTrajectory — Method

DensityTrajectory(system, pulse, initial, goal; algorithm=Tsit5())

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())

Piccolo.Quantum.QuantumTrajectories.KetTrajectory — Type

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.

Piccolo.Quantum.QuantumTrajectories.KetTrajectory — Method

KetTrajectory(system, pulse, initial, goal; algorithm=MagnusGL4())

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: MagnusGL4())

Piccolo.Quantum.QuantumTrajectories.KetTrajectory — Method

KetTrajectory(system, initial, goal, T::Real; drive_name=:u, algorithm=MagnusGL4())

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: MagnusGL4())

Piccolo.Quantum.QuantumTrajectories.MultiKetTrajectory — Type

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.

Piccolo.Quantum.QuantumTrajectories.MultiKetTrajectory — Method

MultiKetTrajectory(system, pulse, initials, goals; weights=..., algorithm=MagnusGL4())

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: MagnusGL4())

Piccolo.Quantum.QuantumTrajectories.MultiKetTrajectory — Method

MultiKetTrajectory(system, initials, goals, T::Real; weights=..., drive_name=:u, algorithm=MagnusGL4())

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: MagnusGL4())

Piccolo.Quantum.QuantumTrajectories.SamplingTrajectory — Type

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

Piccolo.Quantum.QuantumTrajectories.SamplingTrajectory — Method

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)

Piccolo.Quantum.QuantumTrajectories.UnitaryTrajectory — Type

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.

Piccolo.Quantum.QuantumTrajectories.UnitaryTrajectory — Method

UnitaryTrajectory(system, pulse, goal; initial=I, algorithm=MagnusGL4())

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: MagnusGL4())

Piccolo.Quantum.QuantumTrajectories.UnitaryTrajectory — Method

UnitaryTrajectory(system, goal, T::Real; drive_name=:u, algorithm=MagnusGL4())

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: MagnusGL4())

Piccolo.Quantum.Pulses.drive_name — Method

drive_name(traj::AbstractQuantumTrajectory)

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

Piccolo.Quantum.Pulses.duration — Method

duration(traj)

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

Piccolo.Quantum.QuantumTrajectories._add_global_data_to_kwargs — Method

_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

Piccolo.Quantum.QuantumTrajectories._get_control_data — Method

_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

Piccolo.Quantum.QuantumTrajectories._get_control_data — Method

_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)

Piccolo.Quantum.QuantumTrajectories._get_control_data — Method

_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

Piccolo.Quantum.QuantumTrajectories._get_drive_bounds — Method

_get_drive_bounds(sys::QuantumSystem)

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

Piccolo.Quantum.QuantumTrajectories._named_tuple — Method

_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)

Piccolo.Quantum.QuantumTrajectories._sample_times — Method

_sample_times(traj, times::AbstractVector)

Return times as a Float64 vector.

Piccolo.Quantum.QuantumTrajectories._sample_times — Method

_sample_times(traj, N::Int)

Generate N uniformly spaced times for sampling.

Piccolo.Quantum.QuantumTrajectories._sample_times — Method

_sample_times(traj, ::Nothing)

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

Piccolo.Quantum.QuantumTrajectories.extract_pulse — Function

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)

Piccolo.Quantum.QuantumTrajectories.get_goal — Method

get_goal(traj)

Get the goal state/operator from a trajectory.

Piccolo.Quantum.QuantumTrajectories.get_initial — Method

get_initial(traj)

Get the initial state/operator from a trajectory.

Piccolo.Quantum.QuantumTrajectories.get_pulse — Method

get_pulse(traj)

Get the control pulse from a trajectory.

Piccolo.Quantum.QuantumTrajectories.get_solution — Method

get_solution(traj)

Get the ODE solution from a trajectory.

Piccolo.Quantum.QuantumTrajectories.get_system — Method

get_system(traj)

Get the quantum system from a trajectory.

Piccolo.Quantum.QuantumTrajectories.get_systems — Method

get_systems(sampling::SamplingTrajectory)

Get all systems in the sampling trajectory.

Piccolo.Quantum.QuantumTrajectories.get_weights — Method

get_weights(sampling::SamplingTrajectory)

Get the weights for each system.

Piccolo.Quantum.QuantumTrajectories.state_name — Method

state_name(::AbstractQuantumTrajectory)

Get the fixed state variable name for a trajectory type.

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

Piccolo.Quantum.QuantumTrajectories.state_names — Method

state_names(traj::MultiKetTrajectory)

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

Piccolo.Quantum.QuantumTrajectories.state_names — Method

state_names(sampling::SamplingTrajectory)

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

Piccolo.Quantum.QuantumTrajectories.time_name — Method

time_name(::AbstractQuantumTrajectory)

Get the time variable name (always :t).

Piccolo.Quantum.QuantumTrajectories.timestep_name — Method

timestep_name(::AbstractQuantumTrajectory)

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

Piccolo.Quantum.Rollouts._update_system! — Method

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).

Piccolo.Quantum.Rollouts._update_system! — Method

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).

Piccolo.Quantum.Rollouts._update_system! — Method

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).

Piccolo.Quantum.Rollouts._update_system! — Method

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.

Piccolo.Quantum.Rollouts._update_system! — Method

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)

Compute the fidelity between the final unitary and the goal.

rollout!(qtraj::DensityTrajectory, pulse::AbstractPulse; algorithm=Tsit5(), n_points=101)

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_points=101, 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=MagnusGL4(), n_points=101)

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

rollout!(qtraj::KetTrajectory; algorithm=MagnusGL4(), n_points=101, 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=MagnusGL4(), n_points=101)

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

rollout!(qtraj::MultiKetTrajectory; algorithm=MagnusGL4(), n_points=101, 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=MagnusGL4(), n_points=101)

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

rollout!(qtraj::SamplingTrajectory; algorithm=MagnusGL4(), n_points=101, 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=MagnusGL4(), n_points=101)

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: MagnusGL4())
n_points::Int: Number of time points to sample (default: 101)

Example

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

Pulses

Piccolo.Quantum.Pulses.AbstractPulse — Type

AbstractPulse

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

Piccolo.Quantum.Pulses.AbstractSplinePulse — Type

AbstractSplinePulse <: AbstractPulse

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

Piccolo.Quantum.Pulses.CompositePulse — Type

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)]

Piccolo.Quantum.Pulses.CompositePulse — Type

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)

Piccolo.Quantum.Pulses.CubicSplinePulse — Type

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)

Piccolo.Quantum.Pulses.ErfPulse — Type

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

Piccolo.Quantum.Pulses.ErfPulse — Method

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)

Piccolo.Quantum.Pulses.ErfPulse — Method

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)

Piccolo.Quantum.Pulses.GaussianPulse — Type

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

Piccolo.Quantum.Pulses.GaussianPulse — Method

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

Piccolo.Quantum.Pulses.GaussianPulse — Method

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)

Piccolo.Quantum.Pulses.LinearSplinePulse — Type

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)

Piccolo.Quantum.Pulses.LinearSplinePulse — Method

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))

Piccolo.Quantum.Pulses.LinearSplinePulse — Method

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

Piccolo.Quantum.Pulses.LinearSplinePulse — Method

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)

Piccolo.Quantum.Pulses.ZeroOrderPulse — Type

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)

Piccolo.Quantum.Pulses.ZeroOrderPulse — Method

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))

Piccolo.Quantum.Pulses.ZeroOrderPulse — Method

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)

Piccolo.Quantum.Pulses.drive_name — Method

drive_name(pulse::AbstractPulse)

Return the name of the drive variable for this pulse.

Piccolo.Quantum.Pulses.duration — Method

duration(pulse::AbstractPulse)

Return the duration of the pulse.

Piccolo.Quantum.Pulses.get_knot_count — Method

get_knot_count(pulse::AbstractSplinePulse)

Return the number of knots in the spline pulse.

Piccolo.Quantum.Pulses.get_knot_derivatives — Method

get_knot_derivatives(pulse::CubicSplinePulse)

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

Piccolo.Quantum.Pulses.get_knot_times — Method

get_knot_times(pulse::AbstractSplinePulse)

Return the knot times stored in the spline pulse interpolant.

Piccolo.Quantum.Pulses.get_knot_values — Method

get_knot_values(pulse::CubicSplinePulse)

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

Piccolo.Quantum.Pulses.n_drives — Method

n_drives(pulse::AbstractPulse)

Return the number of control drives in the pulse.

Piccolo.Quantum.Pulses.sample — Method

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

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

Piccolo.Quantum.Pulses.sample — Method

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

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

Rollouts

Piccolo.Quantum.Rollouts._update_system! — Function

_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.

Piccolo.Quantum.Rollouts.extract_globals — Function

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.

Piccolo.Quantum.Rollouts.rollout — Function

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

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

Piccolo.Quantum.Rollouts.rollout! — Function

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.

Piccolo.Quantum.Rollouts.unitary_fidelity — Method

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

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

Piccolo.Quantum.Rollouts.update_global_params! — Method

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.

Embedded Operators

Piccolo.Quantum.EmbeddedOperators.AbstractPiccoloOperator — Type

AbstractPiccoloOperator

Union type for operators.

Piccolo.Quantum.EmbeddedOperators.EmbeddedOperator — Type

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.

Piccolo.Quantum.EmbeddedOperators.EmbeddedOperator — Method

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.

Piccolo.Quantum.EmbeddedOperators.EmbeddedOperator — Method

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.

Piccolo.Quantum.EmbeddedOperators.EmbeddedOperator — Method

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

Embed the subspace_operator into a quantum system.

Piccolo.Quantum.EmbeddedOperators.EmbeddedOperator — Method

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.

Piccolo.Quantum.EmbeddedOperators.embed — Method

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

Embed the subspace_operator in the subspace of a larger embedded_operator.

Piccolo.Quantum.EmbeddedOperators.embed — Method

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

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

Piccolo.Quantum.EmbeddedOperators.get_enr_subspace_indices — Method

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.

Piccolo.Quantum.EmbeddedOperators.get_iso_vec_leakage_indices — Function

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.

Piccolo.Quantum.EmbeddedOperators.get_iso_vec_subspace_indices — Function

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.

Piccolo.Quantum.EmbeddedOperators.get_leakage_indices — Function

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.

Piccolo.Quantum.EmbeddedOperators.get_subspace_indices — Function

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.

Piccolo.Quantum.EmbeddedOperators.unembed — Method

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

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

Piccolo.Quantum.EmbeddedOperators.unembed — Method

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

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

Piccolo.Quantum.EmbeddedOperators.unembed — Method

unembed(embedded_op::EmbeddedOperator)

Unembed an embedded operator, returning the original operator.

Lifted Operators

Piccolo.Quantum.LiftedOperators.lift_operator — Function

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 Sums

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

Returns the direct sum of two matrices.

direct_sum(Ã⃗::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.

Isomorphisms

Piccolo.Quantum.Isomorphisms.G — Method

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

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

Piccolo.Quantum.Isomorphisms.H — Method

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

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

Piccolo.Quantum.Isomorphisms.ad_vec — Method

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^*\]

Piccolo.Quantum.Isomorphisms.bloch_to_ket — Method

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

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

Piccolo.Quantum.Isomorphisms.density_to_iso_vec — Method

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

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

Piccolo.Quantum.Isomorphisms.iso — Method

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)\]

Piccolo.Quantum.Isomorphisms.iso_D — Method

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

Returns the isomorphic representation of the Lindblad dissipator L.

Piccolo.Quantum.Isomorphisms.iso_operator_to_iso_vec — Method

iso_operator_to_iso_vec(Ũ::AbstractMatrix{ℝ}) where ℝ <: Real

Convert a real matrix Ũ representing an isomorphism operator into a real vector.

Piccolo.Quantum.Isomorphisms.iso_operator_to_operator — Method

iso_operator_to_operator(Ũ)

Piccolo.Quantum.Isomorphisms.iso_to_ket — Method

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

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

Piccolo.Quantum.Isomorphisms.iso_vec_to_density — Method

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

Returns the density matrix ρ from its isomorphism ρ⃗̃

Piccolo.Quantum.Isomorphisms.iso_vec_to_iso_operator — Method

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

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

Piccolo.Quantum.Isomorphisms.iso_vec_to_operator — Method

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

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

Piccolo.Quantum.Isomorphisms.ket_to_bloch — Method

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

Convert a ket to a Bloch vector representation.

Piccolo.Quantum.Isomorphisms.ket_to_iso — Method

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

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

Piccolo.Quantum.Isomorphisms.mat — Method

mat(x::AbstractVector)

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

Piccolo.Quantum.Isomorphisms.operator_to_iso_operator — Method

operator_to_iso_operator(U)

Piccolo.Quantum.Isomorphisms.operator_to_iso_vec — Method

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

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

Piccolo.Quantum.Isomorphisms.var_G — Method

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.

Quantum Control Problems

Piccolo.Control.QuantumControlProblems.QuantumControlProblem — Type

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)

DirectTrajOpt.Solvers.solve! — Method

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.

Piccolo.Control.QuantumControlProblems.get_trajectory — Method

get_trajectory(qcp::QuantumControlProblem)

Get the NamedTrajectory from the optimization problem.

Piccolo.Control.QuantumControlProblems.sync_trajectory! — Method

sync_trajectory!(qcp::QuantumControlProblem)

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

After optimization, this function:

Extracts the optimized controls from prob.trajectory (unadapting if needed)
Creates a new pulse with those controls via extract_pulse
Re-solves the ODE to get the updated quantum evolution
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

Piccolo.Quantum.Pulses.drive_name — Method

drive_name(qcp::QuantumControlProblem)

Get the control variable name from the quantum trajectory.

Piccolo.Quantum.QuantumTrajectories.get_goal — Method

get_goal(qcp::QuantumControlProblem)

Get the goal state/operator from the quantum trajectory.

Piccolo.Quantum.QuantumTrajectories.get_system — Method

get_system(qcp::QuantumControlProblem)

Get the QuantumSystem from the quantum trajectory.

Piccolo.Quantum.QuantumTrajectories.state_name — Method

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)

Problem Templates

Piccolo.Control.ProblemTemplates.MinimumTimeProblem — Method

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:

Create base problem with SmoothPulseProblem (or similar)
Solve base problem to get feasible solution
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)

Piccolo.Control.ProblemTemplates.SamplingProblem — Method

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

Piccolo.Control.ProblemTemplates.SmoothPulseProblem — Method

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.

Piccolo.Control.ProblemTemplates.SmoothPulseProblem — Method

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.

Piccolo.Control.ProblemTemplates.SmoothPulseProblem — Method

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

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

Piccolo.Control.ProblemTemplates.SplinePulseProblem — Function

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.

Piccolo.Control.ProblemTemplates.SplinePulseProblem — Function

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.

Piccolo.Control.ProblemTemplates.SplinePulseProblem — Method

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

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

Piccolo.Control.ProblemTemplates._ensemble_ket_objective — Method

_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).

Piccolo.Control.ProblemTemplates._final_fidelity_constraint — Method

_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⟩).

Piccolo.Control.ProblemTemplates.add_global_bounds_constraints! — Method

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.

Piccolo.Control.ProblemTemplates.extract_regularization — Method

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).

Piccolo.Control.ProblemTemplates.sampling_state_objective — Method

sampling_state_objective(qtraj, traj, state_sym, Q)

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

Objectives

Piccolo.Control.QuantumObjectives.CoherentKetInfidelityObjective — Method

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)

Piccolo.Control.QuantumObjectives.KetInfidelityObjective — Method

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

Piccolo.Control.QuantumObjectives.KetInfidelityObjective — Method

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

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

Piccolo.Control.QuantumObjectives.LeakageObjective — Method

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.

Piccolo.Control.QuantumObjectives.coherent_ket_fidelity — Method

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

Constraints

Piccolo.Control.QuantumConstraints.FinalCoherentKetFidelityConstraint — Method

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)

Piccolo.Control.QuantumConstraints.LeakageConstraint — Method

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.

Options

Piccolo.Control.Options.PiccoloOptions — Type

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.

Visualizations

Piccolo.Visualizations.QuantumObjectPlots.plot_state_populations — Method

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
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])