This paper presents a class of structure-preserving numerical methods for quantum optimal control problems, based on commutator-free Cayley integrators. Starting from the Krotov framework, we reformulate the forward and backward propagation steps using Cayley-type schemes that preserve unitarity and symmetry at the discrete level. This approach eliminates the need for matrix exponentials and commutators, leading to significant computational savings while maintaining higher-order accuracy. We first recall the standard linear setting and then extend the formulation to nonlinear Schrödinger and Gross-Pitaevskii equations using a Cayley-polynomial interpolation strategy. Numerical experiments on state-transfer problems illustrate that the CF-Cayley method achieves the same accuracy as high-order exponential or Cayley-Magnus schemes at substantially lower cost, especially for longtime or highly oscillatory dynamics. In the nonlinear regime, the structure-preserving properties of the method ensure stability and norm conservation, making it a robust tool for large-scale quantum control simulations. The proposed framework thus bridges geometric integration and optimal control, offering an efficient and reliable alternative to existing exponential-based propagators.
This paper presents a class of structure-preserving numerical methods for quantum optimal control problems, based on commutator-free Cayley integrators. Starting from the Krotov framework, we reformulate the forward and backward propagation steps using Cayley-type schemes that preserve unitarity and symmetry at the discrete level. This approach eliminates the need for matrix exponentials and commutators, leading to significant computational savings while maintaining higher-order accuracy. We first recall the standard linear setting and then extend the formulation to nonlinear Schrödinger and Gross-Pitaevskii equations using a Cayley-polynomial interpolation strategy. Numerical experiments on state-transfer problems illustrate that the CF-Cayley method achieves the same accuracy as high-order exponential or Cayley-Magnus schemes at substantially lower cost, especially for longtime or highly oscillatory dynamics. In the nonlinear regime, the structure-preserving properties of the method ensure stability and norm conservation, making it a robust tool for large-scale quantum control simulations. The proposed framework thus bridges geometric integration and optimal control, offering an efficient and reliable alternative to existing exponential-based propagators.
Optimal control of open quantum systems governed by the Lindblad master equation requires the repeated solution of forward and adjoint evolution equations. For large Hilbert spaces, these computations demand numerical methods that are structure-preserving, accurate, and computationally efficient, while remaining compatible with nonsmooth optimization techniques. In this work, a unified framework for optimal control of Lindblad equations is presented. The continuous optimality system is derived from the Pontryagin maximum principle and solved by a sequential quadratic Hamiltonian (SQH) method. Its numerical realization is based on second-order exponential midpoint propagators for the forward and adjoint Lindblad equations with time-dependent Hamiltonians. Full-rank schemes preserving the Hermitian and positive-semidefinite structure of the density matrix are developed together with low-rank formulations that substantially reduce storage requirements and computational cost. Rigorous error estimates are established for the full- and low-rank forward and adjoint propagators. Numerical experiments confirm the predicted convergence rates, demonstrate the effectiveness of the low-rank approximations, and illustrate the performance of the proposed FREM-SQH and LREM-SQH algorithms for optimal control problems with smooth and nonsmooth control costs.
We consider the numerical propagation of models that combine both quantum and classical degrees of freedom, usually, electrons and nuclei, respectively. We focus, in our computational examples, on the case in which the quantum electrons are modeled with time-dependent density-functional theory, although the methods discussed below can be used with any other level of theory. Often, for these so-called quantum-classical molecular dynamics models, one uses some propagation technique to deal with the quantum part and a different one for the classical equations. While the resulting procedure may, in principle, be consistent, it can however spoil some of the properties of the methods, such as the accuracy order with respect to the time step or the preservation of the geometrical structure of the equations. Few methods have been developed specifically for hybrid quantum-classical models. We propose using the same method for both the quantum and classical particles, in particular, one family of techniques that proves to be very efficient for the propagation of Schrödinger-like equations: the (quasi)-commutator free Magnus expansions. These have been developed, however, for linear systems, yet our problem is nonlinear: formally, the full quantum-classical system can be rewritten as a nonlinear Schrödinger equation, i.e., one in which the Hamiltonian depends on the system itself. The Magnus expansion algorithms for linear systems require the application of the Hamiltonian at intermediate points in a given propagating interval. For nonlinear systems, this poses a problem as this Hamiltonian is unknown due to its dependence on the state. We approximate it by employing a higher order extrapolation using previous steps as input. The resulting technique can then be regarded as a multistep technique or, alternatively, as a predictor corrector formula.
Hamiltonian systems of ordinary and partial differential equations are fundamental mathematical models spanning virtually all physical scales. A critical property for the robustness and stability of computational methods in such systems is the underlying symplectic structure, which preserves geometric properties like phase-space volume over time and energy conservation over an extended period. In this paper, we present quantum algorithms that incorporate symplectic integrators, ensuring the preservation of this key structure. We demonstrate how these algorithms maintain the symplectic properties for both linear and nonlinear Hamiltonian systems. Additionally, we provide a comprehensive theoretical analysis of the computational complexity, showing that our approach offers both accuracy and improved efficiency over classical algorithms. These results highlight the potential application of quantum algorithms for solving large-scale Hamiltonian systems while preserving essential physical properties.
Important nonlinear dynamics, such as those found in plasma and fluid systems, are typically hard to simulate on classical computers. Thus, if fault-tolerant quantum computers could efficiently solve such nonlinear problems, it would be a transformative change for many industries. In a recent breakthrough [Liu et al., PNAS 2021], the first efficient quantum algorithm for solving nonlinear differential equations was constructed, based on a single condition R<1, where R characterizes the ratio of nonlinearity to dissipation. This result, however, is limited to the class of purely dissipative systems with negative log-norm, which excludes application to many important problems. In this work, we correct technical issues with this and other prior analysis, and substantially extend the scope of nonlinear dynamical systems that can be efficiently simulated on a quantum computer in a number of ways. Firstly, we extend the existing results from purely dissipative systems to a much broader class of stable systems, and show that every quadratic Lyapunov function for the linearized system corresponds to an independent R-number criterion for the convergence of the Carleman scheme. Secondly, we extend our stable system results to physically relevant settings where conserved polynomial quantities exist. Finally, we provide extensive results for the class of non-resonant systems. With this, we are able to show that efficient quantum algorithms exist for a much wider class of nonlinear systems than previously known, and prove the BQP-completeness of nonlinear oscillator problems of exponential size. In our analysis, we also obtain several results related to the Poincaré-Dulac theorem and diagonalization of the Carleman matrix, which could be of independent interest.
Abstract Quantum optimal control theory is becoming increasingly crucial as quantum devices become more precise, but the need to quickly optimize these systems classically remains a significant bottleneck in their operation. Here we present a new theoretical quantum control framework for much faster optimization than the state of the art by replacing standard time propagation with a product of short-time propagators, each calculated using the Magnus expansion. The derived formulas for exact series terms and their gradients, based on earlier approximate integrals in a simulation setting, allow us to subsume the high cost of calculating commutators and integrals as an initial overhead. This provides an order of magnitude speedup for quantum control optimization.
We develop a discrete Pontryagin Maximum Principle (PMP) for controlled open quantum systems governed by Lindblad dynamics, and introduce a second--order \emph{contact Lie--group variational integrator} (contact LGVI) that preserves both the CPTP (completely positive and trace--preserving) structure of the Lindblad flow and the contact geometry underlying the discrete PMP. A type--II discrete contact generating function produces a strict discrete contactomorphism under which the state, costate, and cost propagate in exact agreement with the variational structure of the discrete contact PMP.
We apply this framework to the optimal control of a dissipative qubit and compare it with a non--geometric explicit RK2 discretization of the Lindblad equation. Although both schemes have the same formal order, the RK2 method accumulates geometric drift (loss of trace, positivity violations, and breakdown of the discrete contact form) that destabilizes PMP shooting iterations, especially under strong dissipation or long horizons. In contrast, the contact LGVI maintains exact CPTP structure and discrete contact geometry step by step, yielding stable, physically consistent, and geometrically faithful optimal control trajectories.