GitHub

Trajectories

A trajectory bundles a quantum system, a control pulse, and a goal into the central object that problem templates consume.

Overview

Every trajectory type defines:

  • The state$x_k$ and its dynamics $x_{k+1} = \exp(\Delta t_k\, G(\boldsymbol{u}_k))\, x_k$
  • The initial condition$x_1 = x_{\text{init}}$
  • The goal$x_{\text{goal}}$ and associated fidelity metric

Trajectory Types

TypeDynamicsState $x_k$Fidelity $F$$\dim(x_k)$
UnitaryTrajectory$\dot{U} = -iH\,U$$\tilde{U} \in \mathbb{R}^{2d^2}$$\lvert\operatorname{tr}(U_g^\dagger U)\rvert^2 / d^2$$2d^2$
KetTrajectory$\dot{\psi} = -iH\,\psi$$\tilde{\psi} \in \mathbb{R}^{2d}$$\lvert\langle\psi_g\mid\psi\rangle\rvert^2$$2d$
DensityTrajectory$\dot{\rho} = \mathcal{G}\,\text{vec}(\rho)$$\tilde{\rho} \in \mathbb{R}^{d^2}$ (compact)$\operatorname{tr}(\rho_g\,\rho)$$d^2$
MultiKetTrajectory$\dot{\psi}_j = -iH\,\psi_j$multiple $\tilde{\psi}_j$coherent (see below)$2d \times n_{\text{states}}$
MultiDensityTrajectory$\dot{\rho}_j = \mathcal{G}\,\text{vec}(\rho_j)$multiple $\tilde{\rho}_j$weighted average (see below)$d^2 \times M$
SamplingTrajectoryper-system dynamicsper-system statesaverage fidelityvaries

The isomorphic state $x_k$ is always a real vector; see Isomorphisms for the conversions.

UnitaryTrajectory

For synthesizing quantum gates. The propagator satisfies $\dot{U} = G(\boldsymbol{u})\,U$ with $U(0) = I$.

Construction

using Piccolo

# Define system
H_drift = PAULIS[:Z]
H_drives = [PAULIS[:X], PAULIS[:Y]]
sys = QuantumSystem(H_drift, H_drives, [1.0, 1.0])

# Create pulse
T, N = 10.0, 100
times = collect(range(0, T, length = N))
pulse = ZeroOrderPulse(0.1 * randn(2, N), times)

# Create trajectory with goal
U_goal = GATES[:X]
qtraj = UnitaryTrajectory(sys, pulse, U_goal)
UnitaryTrajectory{ZeroOrderPulse{DataInterpolations.ConstantInterpolation{Matrix{Float64}, Vector{Float64}, Vector{Union{}}, Float64}}, SciMLBase.ODESolution{ComplexF64, 3, Vector{Matrix{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, SciMLBase.ODEProblem{Matrix{ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Base.Pairs{Symbol, Vector{Float64}, Nothing, @NamedTuple{tstops::Vector{Float64}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqLinear.MagnusAdapt4, OrdinaryDiffEqCore.InterpolationData{SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Vector{Matrix{ComplexF64}}, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, OrdinaryDiffEqLinear.MagnusAdapt4Cache{Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, Nothing}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}}(QuantumSystem: levels = 2, n_drives = 2, ZeroOrderPulse(Number of drives = 2, T = 10.0), ComplexF64[1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im], ComplexF64[0.0 + 0.0im 1.0 + 0.0im; 1.0 + 0.0im 0.0 + 0.0im], ComplexF64[1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im;;; 0.9948356972119734 - 0.0998277972340243im -0.015546883374002625 + 0.009728353895943075im; 0.015546883374002625 + 0.009728353895943075im 0.9948356972119734 + 0.09982779723402431im;;; 0.9800038616357114 - 0.19863374924825994im -0.007956510413013842 + 0.008588293236771556im; 0.007956510413013838 + 0.00858829323677155im 0.9800038616357112 + 0.19863374924825988im;;; … ;;; -0.9329700936728572 + 0.35794079385866817im -0.031139844811370884 + 0.021806019100029973im; 0.031139844811371217 + 0.02180601910003041im -0.9329700936728608 - 0.35794079385866945im;;; -0.8928039128013926 + 0.4493161670405809im -0.02262148164000127 + 0.022459383136602094im; 0.022621481640001575 + 0.02245938313660257im -0.8928039128013955 - 0.4493161670405825im;;; -0.8432993145417272 + 0.5356907237834225im -0.01853122065217762 + 0.03922127492786102im; 0.018531220652177865 + 0.03922127492786145im -0.8432993145417297 - 0.5356907237834245im])

Solve and Analyze

qcp = SmoothPulseProblem(qtraj, N; Q = 100.0)
solve!(qcp; max_iter = 50)
fidelity(qcp)
0.9997914029235936

Extracting the Pulse

optimized_pulse = get_pulse(qcp.qtraj)
duration(optimized_pulse)
14.262180172049334

KetTrajectory

For state preparation: find controls that map $|\psi_{\text{init}}\rangle \to |\psi_{\text{goal}}\rangle$ (up to global phase). The fidelity is $F = |\langle\psi_{\text{goal}}|\psi(T)\rangle|^2$.

Construction

# Initial and goal states
ψ_init = ComplexF64[1, 0]  # |0⟩
ψ_goal = ComplexF64[0, 1]  # |1⟩

qtraj_ket = KetTrajectory(sys, pulse, ψ_init, ψ_goal)
KetTrajectory{ZeroOrderPulse{DataInterpolations.ConstantInterpolation{Matrix{Float64}, Vector{Float64}, Vector{Union{}}, Float64}}, SciMLBase.ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, Nothing, SciMLBase.ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Base.Pairs{Symbol, Vector{Float64}, Nothing, @NamedTuple{tstops::Vector{Float64}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqLinear.MagnusAdapt4, OrdinaryDiffEqCore.InterpolationData{SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, Nothing, OrdinaryDiffEqLinear.MagnusAdapt4Cache{Vector{ComplexF64}, Vector{ComplexF64}, Matrix{ComplexF64}, Vector{ComplexF64}, Nothing}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}}(QuantumSystem: levels = 2, n_drives = 2, ZeroOrderPulse(Number of drives = 2, T = 10.0), ComplexF64[1.0 + 0.0im, 0.0 + 0.0im], ComplexF64[0.0 + 0.0im, 1.0 + 0.0im], ComplexF64[1.0 + 0.0im 0.9948356972119734 - 0.0998277972340243im … -0.8928039128013926 + 0.4493161670405809im -0.8432993145417272 + 0.5356907237834225im; 0.0 + 0.0im 0.015546883374002625 + 0.009728353895943075im … 0.022621481640001575 + 0.02245938313660257im 0.018531220652177865 + 0.03922127492786145im])

Solve

qcp_ket = SmoothPulseProblem(qtraj_ket, N; Q = 100.0)
solve!(qcp_ket; max_iter = 50)
fidelity(qcp_ket)
0.9999980412691527

MultiKetTrajectory

For gates defined by multiple state mappings with coherent phases. The fidelity enforces phase alignment across all state pairs:

\[F = \left| \frac{1}{n} \sum_{j=1}^{n} \langle \psi_{\text{goal},j} | \psi_j(T) \rangle \right|^2\]

This is strictly harder than per-state fidelity because relative phases must be correct for the mapping to represent a valid gate.

Construction

# Define state pairs: X gate maps |0⟩ → |1⟩ and |1⟩ → |0⟩
ψ0 = ComplexF64[1, 0]
ψ1 = ComplexF64[0, 1]

initial_states = [ψ0, ψ1]
goal_states = [ψ1, ψ0]

qtraj_multi = MultiKetTrajectory(sys, pulse, initial_states, goal_states)
MultiKetTrajectory{ZeroOrderPulse{DataInterpolations.ConstantInterpolation{Matrix{Float64}, Vector{Float64}, Vector{Union{}}, Float64}}, SciMLBase.EnsembleSolution{ComplexF64, 3, Vector{SciMLBase.ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, Nothing, SciMLBase.ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Base.Pairs{Symbol, Vector{Float64}, Nothing, @NamedTuple{tstops::Vector{Float64}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqLinear.MagnusAdapt4, OrdinaryDiffEqCore.InterpolationData{SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, Nothing, OrdinaryDiffEqLinear.MagnusAdapt4Cache{Vector{ComplexF64}, Vector{ComplexF64}, Matrix{ComplexF64}, Vector{ComplexF64}, Nothing}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}}}}(QuantumSystem: levels = 2, n_drives = 2, ZeroOrderPulse(Number of drives = 2, T = 10.0), Vector{ComplexF64}[[1.0 + 0.0im, 0.0 + 0.0im], [0.0 + 0.0im, 1.0 + 0.0im]], Vector{ComplexF64}[[0.0 + 0.0im, 1.0 + 0.0im], [1.0 + 0.0im, 0.0 + 0.0im]], [0.5, 0.5], ComplexF64[1.0 + 0.0im 0.9948356972119734 - 0.0998277972340243im … -0.8928039128013926 + 0.4493161670405809im -0.8432993145417272 + 0.5356907237834225im; 0.0 + 0.0im 0.015546883374002625 + 0.009728353895943075im … 0.022621481640001575 + 0.02245938313660257im 0.018531220652177865 + 0.03922127492786145im;;; 0.0 + 0.0im -0.015546883374002625 + 0.009728353895943075im … -0.022621481640009652 + 0.022459383136611596im -0.018531220652184724 + 0.0392212749278728im; 1.0 + 0.0im 0.9948356972119734 + 0.09982779723402431im … -0.8928039128013945 - 0.4493161670405813im -0.8432993145417292 - 0.5356907237834224im])

Solve

qcp_multi = SmoothPulseProblem(qtraj_multi, N; Q = 100.0)
solve!(qcp_multi; max_iter = 50)
fidelity(qcp_multi)
0.9999654067440789

DensityTrajectory

For open quantum systems governed by the Lindblad master equation. The state $\rho(t)$ evolves as $\dot{\vec\rho} = \mathcal{G}(\boldsymbol{u})\,\vec\rho$ where $\mathcal{G}$ is the Lindbladian superoperator (see Systems).

Internally the state is stored in the compact isomorphism: a real vector of dimension $d^2$ (not $2d^2$) that exploits Hermiticity. The Lindbladian generators are also compacted via $\mathcal{G}_c = P\,\mathcal{G}\,L$ (size $d^2 \times d^2$ instead of $2d^2 \times 2d^2$), giving roughly a 4× speedup in integration. See Isomorphisms for details.

Construction

# Open system with a weak dissipation operator
L = ComplexF64[0.1 0.0; 0.0 0.0]
open_sys = OpenQuantumSystem(
    PAULIS[:Z],
    [PAULIS[:X], PAULIS[:Y]],
    [1.0, 1.0];
    dissipation_operators = [L],
)

# Initial and goal density matrices
ρ_init = ComplexF64[1.0 0.0; 0.0 0.0]  # |0⟩⟨0|
ρ_goal = ComplexF64[0.0 0.0; 0.0 1.0]  # |1⟩⟨1|

T_density, N_density = 10.0, 50
times_density = collect(range(0, T_density, length = N_density))
pulse_density = ZeroOrderPulse(0.1 * randn(2, N_density), times_density)
qtraj_density = DensityTrajectory(open_sys, pulse_density, ρ_init, ρ_goal)
DensityTrajectory{ZeroOrderPulse{DataInterpolations.ConstantInterpolation{Matrix{Float64}, Vector{Float64}, Vector{Union{}}, Float64}}, SciMLBase.ODESolution{ComplexF64, 3, Vector{Matrix{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, SciMLBase.ODEProblem{Matrix{ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Matrix{ComplexF64}, Matrix{ComplexF64}, SciMLBase.NullParameters, Float64}}}, FunctionWrappersWrappers.AllowNonIsBits, FunctionWrappersWrappers.SingleCacheStorage}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Base.Pairs{Symbol, Vector{Float64}, Nothing, @NamedTuple{tstops::Vector{Float64}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), FastBroadcast.Serial}, OrdinaryDiffEqCore.InterpolationData{SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Matrix{ComplexF64}, Matrix{ComplexF64}, SciMLBase.NullParameters, Float64}}}, FunctionWrappersWrappers.AllowNonIsBits, FunctionWrappersWrappers.SingleCacheStorage}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Vector{Matrix{ComplexF64}}, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, OrdinaryDiffEqTsit5.Tsit5Cache{Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), FastBroadcast.Serial}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}}(OpenQuantumSystem: levels = 2, n_drives = 2, ZeroOrderPulse(Number of drives = 2, T = 10.0), ComplexF64[1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im], ComplexF64[0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im], ComplexF64[1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;; 0.999975809887234 + 4.091026692711613e-18im 0.003036117621511072 - 0.0038682612398240196im; 0.0030361176215110727 + 0.0038682612398240183im 2.4190112765953693e-5 - 3.665044666152068e-21im;;; 0.9999042219279081 + 7.51157008156081e-18im 0.005241808967907418 - 0.008260034661078445im; 0.005241808967907418 + 0.008260034661078445im 9.577807209189637e-5 + 4.3353494029362024e-21im;;; … ;;; 0.9850924591104688 + 1.5960545771267064e-16im 0.10714557508726526 + 0.050337099850051324im; 0.10714557508726524 - 0.05033709985005132im 0.014907540889530845 - 1.4858887018517472e-18im;;; 0.9851057675218408 + 1.5949678318386233e-16im 0.11381282286347488 + 0.032151408617280594im; 0.11381282286347487 - 0.03215140861728059im 0.014894232478158917 - 1.6035212447136607e-18im;;; 0.9852884953484585 + 1.5595110127798265e-16im 0.11673253159480565 + 0.013014023654567367im; 0.11673253159480564 - 0.013014023654567355im 0.014711504651540874 - 2.300557702146844e-18im])

Solve and Analyze

The fidelity for density matrices is $F = \operatorname{tr}(\rho_{\text{goal}}\,\rho(T))$.

qcp_density = SmoothPulseProblem(qtraj_density, N_density; Q = 100.0, R = 1e-2)
solve!(qcp_density; max_iter = 150)
fidelity(qcp_density)
0.9948951115506884

MultiDensityTrajectory

For optimizing open quantum systems that must simultaneously steer multiple initial density matrices to their respective targets — e.g. process tomography or channel certification. All density matrices share a single control pulse and Lindblad generator. The fidelity is a weighted average:

\[F = \sum_{j=1}^{M} w_j\, \operatorname{tr}(\rho_{\text{goal},j}\,\rho_j(T)), \qquad \sum_j w_j = 1\]

Internally, M Lindblad ODEs are solved in parallel via DifferentialEquations.EnsembleProblem.

Construction

# Open system: 2-level atom with T₁ decay
L_decay = ComplexF64[0 0; 1 0]  ## |0⟩⟨1| collapse operator (decay rate = 1)
open_sys2 = OpenQuantumSystem(
    PAULIS[:Z],
    [PAULIS[:X], PAULIS[:Y]],
    [1.0, 1.0];
    dissipation_operators = [L_decay],
)

# Prepare two initial → goal pairs (e.g., two rows of a process matrix)
ρ0_a = ComplexF64[1 0; 0 0]  ## |0⟩⟨0|
ρ0_b = ComplexF64[0 0; 0 1]  ## |1⟩⟨1|
ρg_a = ComplexF64[0 0; 0 1]  ## |1⟩⟨1|
ρg_b = ComplexF64[1 0; 0 0]  ## |0⟩⟨0|

T_md, N_md = 10.0, 50
times_md = collect(range(0, T_md, length = N_md))
pulse_md = ZeroOrderPulse(0.1 * randn(2, N_md), times_md)

qtraj_md = MultiDensityTrajectory(
    open_sys2,
    pulse_md,
    [ρ0_a, ρ0_b],
    [ρg_a, ρg_b];
    weights = [0.5, 0.5],
)
MultiDensityTrajectory{ZeroOrderPulse{DataInterpolations.ConstantInterpolation{Matrix{Float64}, Vector{Float64}, Vector{Union{}}, Float64}}, SciMLBase.EnsembleSolution{ComplexF64, 4, Vector{SciMLBase.ODESolution{ComplexF64, 3, Vector{Matrix{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, SciMLBase.ODEProblem{Matrix{ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Matrix{ComplexF64}, Matrix{ComplexF64}, SciMLBase.NullParameters, Float64}}}, FunctionWrappersWrappers.AllowNonIsBits, FunctionWrappersWrappers.SingleCacheStorage}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Base.Pairs{Symbol, Vector{Float64}, Nothing, @NamedTuple{tstops::Vector{Float64}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), FastBroadcast.Serial}, OrdinaryDiffEqCore.InterpolationData{SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Matrix{ComplexF64}, Matrix{ComplexF64}, SciMLBase.NullParameters, Float64}}}, FunctionWrappersWrappers.AllowNonIsBits, FunctionWrappersWrappers.SingleCacheStorage}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Vector{Matrix{ComplexF64}}, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, OrdinaryDiffEqTsit5.Tsit5Cache{Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), FastBroadcast.Serial}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}}}}(OpenQuantumSystem: levels = 2, n_drives = 2, ZeroOrderPulse(Number of drives = 2, T = 10.0), Matrix{ComplexF64}[[1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im], [0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im]], Matrix{ComplexF64}[[0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im], [1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im]], [0.5, 0.5], ComplexF64[1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;; 0.9046027388042209 - 5.276856233320417e-18im 0.007558419772745738 + 0.01214230015820055im; 0.007558419772745738 - 0.012142300158200545im 0.09539726119577917 - 5.951534573722485e-20im;;; 0.8179077432710079 - 7.042916988246711e-18im 0.015380549623257012 + 0.019583050047395877im; 0.015380549623257007 - 0.019583050047395877im 0.182092256728992 - 8.409081043148707e-19im;;; … ;;; 0.00789753868755556 + 1.1538997001793325e-19im -0.0621476826776885 - 0.06203953335200579im; -0.06214768267768847 + 0.06203953335200578im 0.9921024613124426 - 9.18523686107435e-17im;;; 0.005954697255198963 + 3.5690328188389335e-19im -0.06529931749260924 - 0.03922479703813501im; -0.06529931749260921 + 0.039224797038135im 0.9940453027447992 - 9.183854798048359e-17im;;; 0.004513653774228904 + 5.576359177168432e-19im -0.06391094116013823 - 0.01733633567273152im; -0.06391094116013821 + 0.01733633567273152im 0.9954863462257696 - 9.179076848815472e-17im;;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im;;; 0.0002512425946443142 + 1.8403755642349552e-20im -0.00832804389906419 - 0.01348413481808821im; -0.008328043899064195 + 0.013484134818088219im 0.9997487574053556 + 1.4808895076392658e-17im;;; 0.0009467814293497769 - 1.5317404112372242e-20im -0.018632406054771494 - 0.02446798638423689im; -0.018632406054771497 + 0.024467986384236896im 0.9990532185706502 + 3.2065454161475976e-17im;;; … ;;; 0.007805569630956565 + 8.253209842681129e-19im -0.06220901917519297 - 0.06175106430750155im; -0.06220901917519294 + 0.06175106430750157im 0.9921944303690433 + 1.0904596539784408e-16im;;; 0.005875375349078625 + 1.029372647650212e-18im -0.06530122496900076 - 0.038943086036403im; -0.06530122496900073 + 0.03894308603640302im 0.9941246246509212 + 1.0786107614693321e-16im;;; 0.004446010849747049 + 1.2304937460097876e-18im -0.0638588456654172 - 0.017072331273075062im; -0.06385884566541719 + 0.017072331273075097im 0.9955539891502531 + 1.0849637739466801e-16im])

Rollout and Fidelity

qtraj_md_out = rollout(qtraj_md)
fidelity(qtraj_md_out)
0.4999661785377583
Note

Built-in optimization support for MultiDensityTrajectory (via SmoothPulseProblem) requires a Piccolissimo integrator that provides per-density Lindbladian dynamics. Rollout and fidelity evaluation work out of the box with any pulse type.

When to Use vs DensityTrajectory

ScenarioUse
Single initial stateDensityTrajectory
Process tomography / multiple initial statesMultiDensityTrajectory
Weighted importance across initial statesMultiDensityTrajectory with custom weights

SamplingTrajectory

For robust optimization over parameter variations. Created internally by SamplingProblem, this trajectory type stores multiple system variants sharing a single control pulse. The objective averages fidelity across all samples:

\[\bar{\ell} = \frac{1}{S}\sum_{s=1}^{S} \ell\!\bigl(x_N^{(s)},\, x_{\text{goal}}\bigr)\]

systems = [sys_nominal, sys_high, sys_low]
qcp_robust = SamplingProblem(qcp_base, systems)

Common Operations

Extracting and Saving Pulses

After solving, extract the optimized pulse with get_pulse and persist it with save:

optimized_pulse = get_pulse(qcp.qtraj)
duration(optimized_pulse)
14.262180172049334

Save and reload using the pulse-aware JLD2 helpers:

save("my_gate.jld2", optimized_pulse)
saved_pulse = load_pulse("my_gate.jld2")
qtraj_warm = UnitaryTrajectory(sys, saved_pulse, U_goal)
qcp_warm = SmoothPulseProblem(qtraj_warm, N; Q=1000.0)
solve!(qcp_warm; max_iter=50)  # converges faster from a good initial guess

See Saving and Loading Pulses for the full guide.

Named Trajectory Integration

After solving, the NamedTrajectory stores the NLP decision vector at each of the $N$ knot points:

traj = get_trajectory(qcp)

# Controls at timestep k
u_1 = traj[1][:u]
u_1

# All timesteps
Δts = get_timesteps(traj)
length(Δts)
100

Internal Representation

States in the trajectory are real isomorphic vectors. Convert back to complex form with:

U = iso_vec_to_operator(traj[:Ũ⃗][:, end])     # unitary
ψ = iso_to_ket(traj[:ψ̃][:, end])               # ket
ρ = compact_iso_to_density(traj[:ρ⃗̃][:, end])   # density matrix

See Isomorphisms for details.

Best Practices

1. Match Pulse Type to Problem

# For SmoothPulseProblem
pulse = ZeroOrderPulse(controls, times)
qtraj = UnitaryTrajectory(sys, pulse, U_goal)
qcp = SmoothPulseProblem(qtraj, N)  # ✓

# For SplinePulseProblem
pulse = CubicSplinePulse(controls, tangents, times)
qtraj = UnitaryTrajectory(sys, pulse, U_goal)
qcp = SplinePulseProblem(qtraj)  # ✓

2. Initialize with Reasonable Controls

# Scale by drive bounds
max_amp = 0.1 * 1.0
initial_controls = max_amp * randn(2, N)
extrema(initial_controls)
(-0.24976391118628466, 0.2827190729912068)

See Also

This page was generated using Literate.jl.