Finite elements and moving asymptotes accelerate quantum optimal control -- FEMMA

Quantum optimal control is central to designing spin manipulation pulses. Gradient-based pulse optimization can be facilitated by either accelerating gradient evaluation or enhancing the convergence rate. In this work, we accelerated single-spin optimal control by combining the finite element method with the method of moving asymptotes. By treating discretized time as spatial coordinates, the Liouville - von Neumann equation was reformulated as a linear system, efficiently yielding a joint solution of the spin trajectory and control gradient. The method of moving asymptotes, relying on the ensemble fidelities and gradients, achieves rapid convergence for a target fidelity of 0.995.

💡 Research Summary

The paper introduces a novel framework for accelerating quantum optimal control of single‑spin magnetic‑resonance pulses by merging the finite element method (FEM) with the method of moving asymptotes (MMA). The authors reinterpret the discretized time axis as a spatial domain and apply a Galerkin‑type FEM to the Liouville–von Neumann (LvN) equation, converting it into a linear system K·α = f. Here, α contains the density‑matrix values at each temporal node, K is a complex‑valued stiffness matrix assembled from element contributions, and f encodes the initial state. By using linear, quadratic, or Hermite shape functions, the method yields a banded sparse matrix that can be solved efficiently with standard sparse solvers.

A key advantage of this formulation is that the gradient of the fidelity with respect to the control amplitudes can be obtained with essentially the same computational effort as the forward propagation. The gradient expression involves K⁻¹·(∂f/∂x – (∂K/∂x)·α). Because the control waveform is piecewise‑constant, ∂K/∂x is non‑zero only on the element where the control acts, allowing a local evaluation. The adjoint vector λ is defined by the linear system Kᵀλ = (∂η/∂α)ᵀ, and once λ is computed the full gradient follows from λᵀ·(∂f/∂x – (∂K/∂x)·α). Consequently, the cost of gradient evaluation scales linearly with the number of time steps and is comparable to a single LvN solve, dramatically reducing the overhead typical of GRAPE‑type algorithms that require separate forward and backward propagations.

To enforce realistic hardware constraints (e.g., maximum RF amplitude) and to promote smooth pulse shapes, the authors embed a Helmholtz filter into the workflow. The filter is itself a linear FEM problem K_h·x_s = f_h, where x_s denotes the smoothed control variables. After filtering, a hyperbolic‑tangent scaling x = (1 – e^{–κx_s})/(1 + e^{–κx_s}) maps the smoothed amplitudes into the admissible interval