Time-Dependent Quantum Control Gradient

• Time-dependent quantum control gradients are derivatives of optimization objectives with respect to dynamic controls, enabling precise tuning of quantum systems.
• They employ a forward propagation and backward adjoint structure to compute gradients accurately across closed, open, and non-Hermitian quantum dynamics.
• Advanced computational strategies, including exact slice-wise derivatives and hybrid automatic differentiation, enhance efficiency and accuracy in quantum control simulations.

Searching arXiv for recent and foundational papers on time-dependent quantum control gradients, GRAPE, adjoints, and open-system control. A time-dependent quantum control gradient is the derivative of a control objective with respect to time-dependent control variables in a dynamical quantum model. In the literature covered here, those variables include coherent amplitudes in Hamiltonians of the form $H(t)=H_0+\sum_k u_k(t)H_k$, incoherent controls entering decoherence rates $\gamma_k(t)$, basis coefficients in function-space parametrizations, and parameters of neural or spline-based control ansätze. The central technical problem is to differentiate a time-ordered propagation—unitary, non-Hermitian, or Lindbladian—while preserving the numerical structure of the forward solver. Across GRAPE, Krotov, discrete adjoints, semi-automatic differentiation, and projected-gradient formulations, the recurring pattern is a forward propagation, a backward adjoint or costate propagation, and a local contraction between the two through an exact or approximate derivative of a single time-step propagator (Chowdhury et al., 2023).

1. Dynamical setting and objective functionals

Time-dependent quantum control gradients are defined relative to a controlled evolution law and an optimization functional. For closed systems, the standard model is the time-dependent Schrödinger equation with

$H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$

or, in operator form,

$$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$

with the solution expressed as a time-ordered exponential (Lucarelli, 2016). In open systems, the corresponding controlled dynamics are written as time-dependent Lindblad or GKSL equations,

$$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$

or, with time-dependent decoherence rates,

$$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$

For radical pairs, spin-selective recombination is represented in Haberkorn form and embedded into a non-Hermitian generator $A(t)=H(t)-iK$, so that

$$\dot{\rho}(t)=-i\big(A(t)\rho(t)-\rho(t)A^\dagger(t)\big)$$

and $\rho(t)=U(t,0)\rho(0)U^\dagger(t,0)$ with a time-ordered non-Hermitian propagator (Chowdhury et al., 2023).

The objective may be a terminal fidelity, a gate figure of merit, a time-integrated observable, or a multi-term Bolza functional. Representative examples include the transition or gate objectives

$$J(u)=|\langle \psi_{\rm target}|\psi(T)\rangle|^2,\qquad J(u)=\operatorname{Tr}[O\,\rho(T)],$$

the radical-pair reaction yields

$\gamma_k(t)$0

and the TDKS objective

$\gamma_k(t)$1

(Chowdhury et al., 2023, Sprengel et al., 2017). Open-system gate synthesis with coherent and incoherent controls uses a Hilbert–Schmidt discrepancy over four basis states,

$\gamma_k(t)$2

while function-space approaches optimize a phase-invariant fidelity

$\gamma_k(t)$3

(Petruhanov et al., 2023, Lucarelli, 2016).

A common misconception is that the gradient is defined only for terminal-time fidelity. The cited literature treats terminal objectives, running costs, time-integrated yields, state-tracking objectives, and non-analytic gate measures within a unified differentiation framework (Chowdhury et al., 2023, Goerz et al., 2022).

2. Adjoint-state structure and continuous-time gradient formulas

The standard gradient structure is an adjoint or costate formulation. For open Lindbladian dynamics,

$\gamma_k(t)$4

the sensitivity $\gamma_k(t)$5 obeys

$\gamma_k(t)$6

and the objective derivative is $\gamma_k(t)$7. Introducing the co-state $\gamma_k(t)$8 by

$\gamma_k(t)$9

yields the continuous-time adjoint formula

$H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$0

For control-linear Hamiltonians and control-dependent dissipators,

$H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$1

(He et al., 2024).

Closed-system adjoint formulas have the same structure. In a Schrödinger model with terminal expectation or fidelity, the co-state evolves backward under the same Hamiltonian, and the continuous-time gradient is

$H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$2

For control coefficients in a basis expansion $H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$3, the parameter gradient is

$H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$4

(Li et al., 2023).

TDKS optimal control introduces a distinct adjoint because the Hamiltonian depends on the density through $H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$5 and $H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$6. The reduced gradient is

$H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$7

where $H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$8 is the $H(t)=H_0+\sum_{k=1}^K u_k(t)H_k,$9-Riesz representative of $$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$0, equivalently

$$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$1

with the stated boundary conditions when working in $$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$2 (Sprengel et al., 2017).

These formulations show that the gradient is not an abstract derivative of a scalar cost alone; it is a bilinear contraction between a forward-propagated quantum object and a backward-propagated dual object. This remains true for unitary dynamics, Lindbladians, non-Hermitian radical-pair models, and density-dependent TDKS equations.

3. Time discretization, propagator derivatives, and exact slice-wise formulas

In practice, the control is discretized. With piecewise-constant controls on a grid $$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$3 split into $$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$4 steps of size $$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$5, one defines step propagators such as

$$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$6

or, for radical pairs in Liouville space,

$$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$7

The discrete adjoint recursion then yields a slice-wise gradient. In Liouville space,

$$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$8

with backward recursion

$$\dot{U}(t)=-\frac{i}{\hbar}\left(H_d+\sum_{j=1}^M \Omega_j(t)H_j\right)U(t),$$9

(Chowdhury et al., 2023).

The central technical object is the Fréchet derivative of a matrix exponential. For a perturbation $$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$0 of a generator $$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$1,

$$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$2

Hence, for a slice propagator,

$$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$3

The same Wilcox or Fréchet formula appears in open-system GRAPE with coherent and incoherent controls, in two-level GKSL gate generation, and in function-space variants when differentiating step exponentials (Petruhanov et al., 2023, Petruhanov et al., 2023).

Several papers emphasize that the common first-order approximation

$$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$4

is not the exact slice derivative. Exact discrete adjoint methods instead derive gradients of the discrete objective compatible with the numerical integrator. For Störmer–Verlet discretization, the discrete gradient is computed exactly at the cost of two Schrödinger solves, independently of the number of control parameters (Petersson et al., 2020). In Crank–Nicolson form, the discrete gradient for a scalar control is written as

$$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$5

(Riahi et al., 2016).

Radical-pair control extends this exact slice logic to non-Hermitian evolution. The paper’s “central single time-step propagator” evaluates each gradient contribution using the forward state at the beginning of the slice, the backward costate at the end of the slice, and the exact Fréchet derivative of the central propagator. The adjoint must be taken with respect to the Hilbert–Schmidt inner product, not by naive Hermitian conjugation of a non-Hermitian Hamiltonian (Chowdhury et al., 2023).

4. Computational realizations and acceleration strategies

A large part of the modern literature is devoted to making the gradient computation tractable. One strategy is structural factorization of the propagator derivative. For generators split as $$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$6, symmetric Trotter–Suzuki splittings approximate the step propagator: $$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$7 for second-order Strang splitting, and

$$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$8

for fourth-order Yoshida construction, with the Fréchet derivative propagated through the product rule (Chowdhury et al., 2023). In TDKS, forward and adjoint equations are discretized by Strang splitting, giving second-order convergence in time and spectral convergence in space for smooth potentials (Sprengel et al., 2017).

A second strategy is to reduce the control dimensionality. In GRAFS, each control is expanded in Slepian sequences,

$$\dot{\rho}(t)=-i[H(t),\rho(t)]+\sum_j \mathcal{D}[L_j]\rho(t),$$9

and the coefficient-space gradient is the projected time-grid gradient

$$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$0

This parameterization restricts the search space to band-limited controls and recovers the time-bandwidth quantum speed-limit scaling (Lucarelli, 2016). Related dimensionality reduction appears in B-spline-plus-carrier parameterizations, where the number of control parameters is independent of the number of integration steps and the discrete adjoint remains exact (Petersson et al., 2020).

A third strategy is to reduce memory or sampling cost. Quantum trajectories and automatic differentiation replace density-matrix propagation by Monte Carlo wave-function evolution. The gradient is estimated as

$$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$1

with reverse-mode AD through no-jump and jump updates. An improved-sampling algorithm exploits rare jumps by computing the no-jump gradient separately and combining it with a weighted average over jump trajectories (Abdelhafez et al., 2019). Semi-automatic differentiation takes a different route: it applies AD only to the terminal functional, while propagation derivatives are evaluated analytically via a modified GRAPE scheme and a block “gradient generator,” avoiding the prohibitive memory and runtime overhead of full AD through the entire propagation (Goerz et al., 2022).

Time parallelization provides a fourth acceleration. The Intermediate State Method decomposes the horizon into subintervals, defines intermediate states

$$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$2

and proves a segment-wise gradient decomposition

$$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$3

This enables parallel evaluation of local gradient contributions and yields computational time approximately divided by the number of available processors when based on a gradient method (Riahi et al., 2016).

5. Methodological variants across closed, open, constrained, and learned control

The time-dependent quantum control gradient is not tied to a single optimization method. In GRAPE, the gradient is used in direct ascent on a time grid or in a reduced basis. In Krotov’s method, the same forward–backward structure yields an on-the-fly update,

$$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$4

with a second-order correction required for non-convex final-time functionals such as the geometric perfect-entangler objective (Goerz et al., 2015).

For constrained controls, the Gradient Projection Method formulates the gradient directly in terms of quantum objects and enforces pointwise bounds exactly by projection. For closed systems,

$$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$5

and for open systems with simultaneous coherent and incoherent controls,

$$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$6

Heavy-ball two-step variants add inertial terms while preserving exact satisfaction of bounds (Morzhin et al., 2024).

Differentiable programming and semi-automatic differentiation generalize the notion of a control variable beyond piecewise-constant amplitudes. In one approach, a neural control agent outputs $$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$7 from the current state and previous controls, and gradients are obtained by automatic differentiation through both the neural network and the ODE solver. In continuous-time notation, the chain rule takes the form

$$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$8

while the functional derivative $$\dot{\rho}(t)= -\,i\,[H(t),\,\rho(t)] + \sum_{k} \gamma_k(t)\Big(L_k \rho(t) L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k,\rho(t)\}\Big).$$9 remains the familiar adjoint contraction (Schäfer et al., 2020).

A related point concerns non-analytic objectives. Semi-automatic differentiation explicitly separates time propagation from terminal functional evaluation so that objectives involving eigenvalue decompositions, complex logarithms, and max operations—such as gate concurrence—can still be optimized with gradient-based methods. This does not remove non-differentiability at kinks, but it allows branch-dependent subgradients to be used effectively in practice (Goerz et al., 2022).

6. Applications, comparisons, and open issues

Applications of time-dependent quantum control gradients span radical-pair chemistry, TDDFT, gate synthesis, molecular orientation, Bose–Einstein condensates, and quantum readout. Radical-pair control demonstrates time-global optimization of spin-selective recombination yields in the low magnetic field regime, including toy models, a realistic exciplex-forming donor–acceptor system comprising 16 nuclear spins, and spin-biological examples (Chowdhury et al., 2023). TDKS control provides existence and characterization results for bilinear control of multi-electron systems within TDDFT, together with a nonlinear conjugate gradient implementation (Sprengel et al., 2017). Open-system GRAPE with coherent and incoherent controls treats time-dependent decoherence rates and even analytical single-qubit Bloch-space gradients obtained through Cardano diagonalization of a $A(t)=H(t)-iK$0 matrix (Petruhanov et al., 2023). Trajectory-based methods demonstrate state transfer, gate-time reduction, and readout optimization while often reducing the computational cost relative to density-matrix optimization (Abdelhafez et al., 2019).

The literature also draws several careful comparisons. One is between time-local and time-global optimization. For radical pairs, time-local optimization can be effective in high bias fields or for single-time targets, but low-field reaction yields are time-integrated quantities, and controls that temporarily decrease instantaneous singlet probability may still improve the net yield. In that setting, time-global GRAPE outperforms time-local optimization, especially for asymmetric recombination and non-Hermitian dynamics (Chowdhury et al., 2023). Another comparison concerns continuous control spaces versus time-grid parameterizations: TDKS optimization is posed in $A(t)=H(t)-iK$1, with the control gradient expressed through an $A(t)=H(t)-iK$2-Riesz representative rather than only an $A(t)=H(t)-iK$3 pointwise derivative (Sprengel et al., 2017). A third comparison concerns full AD versus hybrid analytic differentiation: semi-automatic differentiation and discrete adjoints are motivated by the observation that naive full AD through long quantum propagations incurs prohibitive memory and runtime overhead (Goerz et al., 2022, Petersson et al., 2020).

Two broader directions extend the subject beyond classical simulation. One line develops hybrid quantum-classical algorithms for open quantum systems, where a quantum procedure simulates time-dependent Lindblad evolution and provides access to gradients via quantum gradient estimation, while a classical procedure updates the controls (He et al., 2024). Another gives efficient quantum algorithms for optimal control under time-dependent Schrödinger dynamics, using sparse-oracle Hamiltonian simulation and fast gradient estimation to obtain exponential speedup in the number of qubits relative to classical algorithms that scale in $A(t)=H(t)-iK$4 (Li et al., 2023). These works preserve the same underlying gradient logic—state propagation, adjoint structure, and parameter differentiation—but relocate the computational bottleneck.

A persistent misconception is that automatic differentiation alone solves the gradient problem. The cited work suggests a more qualified picture. Automatic differentiation can compute pathwise derivatives through trajectories, neural policies, or non-analytic terminal functionals, but accurate and efficient optimization still depends on the propagation scheme, the correct adjoint, control parameterization, and constraint handling. Another misconception is that open-system gradients are a minor modification of unitary GRAPE. In fact, control-dependent dissipators, non-Hermitian generators, density-dependent mean fields, and positivity constraints can materially change both the derivation and the numerical realization of the gradient (Petruhanov et al., 2023, Chowdhury et al., 2023, Sprengel et al., 2017).

In this sense, the time-dependent quantum control gradient is best understood not as a single formula but as a family of structurally related derivatives: all are built from forward propagation, backward dual propagation, and exact or controlled approximation of the derivative of the time-step evolution map, but each dynamical setting fixes its own adjoint, inner product, constraint set, and computational realization.