This work introduces the High-Order Hermite Optimization (HOHO) method, an open-loop discrete adjoint method for quantum optimal control. Our method is the first of its kind to efficiently compute exact (discrete) gradients when using continuous, parameterized control pulses while solving the forward equations (e.g. Schrodinger’s equation or the Linblad master equation) with an arbitrarily high-order Hermite Runge-Kutta method. The HOHO method is implemented in QuantumGateDesign.jl, an open-source software package for the Julia programming language, which we use to perform numerical experiments comparing the method to Juqbox.jl. For realistic model problems we observe speedups up to 775x.

It is well known that numerical solutions of Helmholtz boundary-value problems suffer from pollution (dispersion) errors that require much finer grids than expected for high frequencies. In this article a dispersion error analysis is presented that leads to an explicit rule-of-thumb formula to estimate the points-per-wavelength required for a p-th order accurate finite difference approximation for Helmholtz problems that accounts for pollution errors. The rule of thumb only depends on the wave-number (or frequency) of the periodic forcing, the size of the domain, and the desired relative error tolerance. Several numerical examples on simple and complex domains and with closed or open boundaries show that the rule of thumb is a useful predictor of the error. Having such a rule is important in practice to either avoid under-resolving the problem and obtaining the wrong answer, or over-resolving the problem resulting in a much more expensive computation.

Fast and accurate classical simulation of quantum systems is a central challenge in the design and control of quantum computers, but the highly oscillatory dynamics of these systems severely limit the efficiency of standard numerical methods. To address this, we adapt Filon quadrature for oscillatory integrals into two numerical methods, called Filon and Controlled Filon, for solving linear systems of ODEs with highly oscillatory solutions. We tailor both methods for efficient implementation in controlled quantum systems, and the Controlled Filon method additionally accounts for the oscillatory structure of the control pulses. We show by numerical experiments that these methods significantly reduce the computational cost of accurately simulating systems of superconducting transmon qubits by decreasing the number of timesteps needed to reach a given level of precision, with only a modest increase in the cost per timestep. For a realistic simulation of the dynamics of a CNOT gate, the Controlled Filon method is the most efficient method tested at every target accuracy, outperforming the best Hermite method by up to 6x and the Hermite method of the same order by up to 500x.

We consider the numerical solution of the wave equation in materials with rapidly varying coefficients, and time harmonic sources. For these problems, direct discretization is prohibitively costly, and instead multiscale methods are used. There are several multiscale methods that directly discretize in the frequency domain. In this work we instead start in the time-domain and combine a finite difference Heterogeneous Multiscale Method (HMM) for the wave equation with the WaveHoltz method. Each WaveHoltz iteration marches the wave equation towards the time-periodic Helmholtz solution. The advantages of the WaveHoltz method relative to traditional Helmholtz solvers carry over directly to the multiscale problems considered here. Since, in addition, the time-domain solver does not artificially impose boundary conditions on the micro-scale problems, no boundary errors from the micro-scale problems are present in the homogenized frequency domain solution.

We propose a low-rank method for solving the Helmholtz equation. Our approach is based on the WaveHoltz method, which computes Helmholtz solutions by applying a time-domain filter to the solution of a related wave equation. The wave equation is discretized by high-order multiblock summation-by-parts finite differences. In two dimensions we seek to compress the solution in matrix form, and in three dimensions using tensor trains. To control rank growth we use step-truncation during time stepping and a low-rank Anderson acceleration for the WaveHoltz fixed point iteration. We have carried out extensive numerical experiments demonstrating the convergence and efficacy of the iterative scheme for free- and half-space problems in two and three dimensions with constant and piecewise constant wave speeds.

Solving the Helmholtz equation with iterative methods is challenging because of its indefinite nature and highly oscillatory solutions. The WaveHoltz algorithm mitigates the difficulties of solving the Helmholtz equation by repeatedly solving the wave equation over short periods in time. In this paper, we prove that for stable semi-discretizations of the wave equation, the WaveHoltz iteration converges to an approximate solution of the corresponding frequency-domain problem, provided one exists. We present numerical examples in one and two dimensions using finite difference and discontinuous Galerkin discretizations illustrating these convergence results.

Numerical simulation of quantum computing hardware and open quantum systems governed by the Lindblad equation is challenging due to the high dimensionality of the density matrix and the need to preserve fundamental physical properties. In our previous work, we developed an arbitrary-order, low-rank, completely positive and trace preserving (CPTP) method for the Lindblad equation with time-dependent Hamiltonians by nested Picard iteration (NPI). In this work, we develop Gregory NPI schemes, which are CPTP schemes constructed by Gregory-type quadrature on equispaced nodes. The methods, which are of order up to nine, substantially reduce the computational cost compared to our previously proposed NPI schemes with Gaussian quadrature rules, while retaining high-order accuracy and structure preservation. We analyze the stability of the resulting scheme for a physics-based test equation. Numerical experiments verify the convergence of the method and demonstrate the effectiveness of the low-rank approximation. We study the performance of a previously constructed CNOT gate for both closed and open quantum systems.

We propose a family of low-rank, completely positive and trace preserving schemes for the Lindblad equation, a common model for open quantum systems. Low-rank representation is employed at two levels: the density matrix is factorized into the product of tall-skinny matrices, and the columns of these matrices are further represented using the tensor train (TT) format, also know as matrix product states (MPS). This two-level low-rank format fits naturally into our existing Kraus is King scheme (arXiv:2409.08898v2 [math.NA]) for the Lindblad equation, whose underlying operations are arithmetic on the columns of the tall-skinny matrices. We show how these operations can be performed efficiently in the TT/MPS format, with particular emphasis on density matrix rank-truncation. We conclude with extensive numerical experiments demonstrating the convergence of this scheme and its efficiency in simulating systems with up to 10^19 degrees of freedom using only modest compute resources.

In this paper, we develop a framework for designing arbitrary high order low-rank schemes for the Lindblad equation with time-dependent Hamiltonians. Our approach is based on nested Picard iterative integrators (NPI) and results in schemes in Kraus form that are completely positive and trace preserving (CPTP). The schemes are amenable to low rank formulations, making them suitable for problems where the matrix rank of the density matrix is small.

WAVES 2026, 17th International Conference on Mathematical and Numerical Aspects of Wave Propagation Concordia University (John Molson Building), Montreal, Canada, June 22-26, 2026 WAVES 2026 is the seventeenth meeting in a long-running biennial series that has, throughout its history, alternated between Europe and North America to advance the mathematical and numerical study of wave propagation. Hosted at the John Molson Building of Concordia University in Montreal from June 22 to 26, 2026, the conference continues this tradition as a leading international forum where theory, computation, and application meet. The scientific program spans the full breadth of mathematical and numerical techniques for wave phenomena, from the modeling and analysis of the governing partial differential equations to the design, analysis, and implementation of efficient computational methods. Representative themes include acoustic, electromagnetic, elastic, and seismic wave propagation; scattering and inverse problems; high-frequency and asymptotic methods; finite element, boundary integral, and time-domain discretizations; and absorbing boundary conditions and domain truncation, with applications reaching across the physical sciences and engineering.

We present a numerical study of eigenvector deflation as a means of accelerating the WaveHoltz method for solving the Helmholtz equation. For energy-conserving (Dirichlet or Neumann) boundary conditions the WaveHoltz fixed-point iteration converges slowly at high frequency, requiring approximately 𝒪(ω^2d) iterations in d dimensions. We show that deflating the eigenvectors whose eigenvalues lie nearest the driving frequency substantially reduces iteration counts, and we examine two ways of incorporating the eigenvectors: direct eigenvector deflation (DEVD), in which the forcing and iterate are projected against the deflation set, and augmented-Krylov eigenvector deflation (AUKED) using deflated conjugate gradient (DCG), augmented GMRES (AGMRES), and augmented (recycled) BICGSTAB (ABICGSTAB). The required eigenpairs can be computed efficiently with the EigenWave approach, and we demonstrate, in two dimensions, that when the number of deflation vectors grows quadratically with ω the asymptotic convergence rate remains essentially constant. Because the eigenvectors on structured grids are naturally represented as matrices, we further apply SVD-based compression to reduce their storage. Numerical experiments on single curvilinear grids discretized with summation-by-parts operators, and on overset grids illustrate the robustness and efficiency of the approach, with the deflated solver breaking even against the undeflated solver after as few as two right-hand sides, when accounting for the cost of precomputing the eigenvectors.

An algorithm, named EigenWave, is described to compute eigenvalues and eigenvectors of elliptic boundary value problems. EigenWave is based on the recently developed WaveHoltz scheme and solves a related time-dependent wave equation as part of an iteration involving an initial condition of the time-dependent problem. At each iteration, the solution to the wave equation is filtered in time resulting in updated initial data. The time-filter is designed to promote the contributions of eigenmodes in the solution near a chosen target frequency (target eigenvalue). The ability to choose an arbitrary target frequency enables the computation of eigenvalues anywhere in the spectrum, without the need to invert an indefinite matrix, as is common with other approaches. Furthermore, the iteration can be embedded within a matrix-free Arnoldi algorithm, which enables the efficient computation of multiple eigenpairs near the target frequency. For efficiency, the time-dependent wave equation can be solved with implicit time-stepping and only about 10 time-steps per-period are needed, independent of the mesh spacing. When the (definite) implicit time-stepping equations are solved with a multigrid algorithm, the cost of the resulting EigenWave scheme scales linearly with the number of spatial grid points N as the mesh is refined, giving an optimal O(N) algorithm. The approach is demonstrated by finding eigenpairs of the Laplacian in complex geometry using finite-difference approximations on overset grids. Results in two and three space dimensions are presented using second-order and fourth-order accurate approximations.

We develop efficient and high-order accurate solvers for the Helmholtz equation on complex geometry. The solvers are based on the WaveHoltz algorithm, which computes solutions of the Helmholtz equation by time-filtering solutions of the wave equation. The approach avoids the need to invert an indefinite matrix which can cause convergence difficulties for many iterative solvers for indefinite Helmholtz problems. Complex geometry is treated with overset grids which use Cartesian grids throughout most of the domain together with curvilinear grids near boundaries. The basic WaveHoltz fixed-point iteration is accelerated using GMRES and also by a deflation technique using a set of precomputed eigenmo des. The solution of the wave equation is solved efficiently with implicit time-stepping using as few as five time-steps per period, independent of the mesh size. When multigrid is used to solve the implicit time-stepping equations, the cost of the resulting WaveHoltz scheme scales linearly with the total number of grid points N (at fixed frequency) and is thus optimal in CPU time and memory usage as the mesh is refined. Numerical results are given for problems in two and three dimensions, to second- and fourth-order accuracy, and they show the potential of the approach to solve a wide range of large-scale problems.

This work proposes a new class of preconditioners for the low rank Generalized Minimal Residual Method (GMRES) for multiterm matrix equations arising from implicit timestepping of linear matrix differential equations. We are interested in computing low rank solutions to matrix equations, e.g. arising from spatial discretization of stiff partial differential equations (PDEs). The low rank GMRES method is a particular class of Krylov subspace method where the iteration is performed on the low rank factors of the solution. Such methods can exploit the low rank property of the solution to save on computational and storage cost. Of critical importance for the efficiency and applicability of the low rank GMRES method is the availability of an effective low rank preconditioner that operates directly on the low rank factors of the solution and that can limit the iteration count and the maximal Krylov rank. The preconditioner we propose here is based on the basis update and Galerkin (BUG) method, resulting from the dynamic low rank approximation. It is a nonlinear preconditioner for the low rank GMRES scheme that naturally operates on the low rank factors. Extensive numerical tests show that this new preconditioner is highly efficient in limiting iteration count and maximal Krylov rank. We show that the preconditioner performs well for general diffusion equations including highly challenging problems, e.g. high contrast, anisotropic equations. Further, it compares favorably with the state of the art exponential sum preconditioner. We also propose a hybrid BUG - exponential sum preconditioner based on alternating between the two preconditioners.

We propose a novel Hermite-Taylor correction function method to handle embedded boundary and interface conditions for Maxwell's equations. The Hermite-Taylor method evolves the electromagnetic fields and their derivatives through order m in each Cartesian coordinate. This makes the development of a systematic approach to enforce boundary and interface conditions difficult. Here we use the correction function method to update the numerical solution where the Hermite-Taylor method cannot be applied directly. Time derivatives of boundary and interface conditions, converted into spatial derivatives, are enforced to obtain a stable method and relax the time-step size restriction of the Hermite-Taylor correction function method. The proposed high-order method offers a flexible systematic approach to handle embedded boundary and interface problems, including problems with discontinuous solutions at the interface. This method is also easily adaptable to other first order hyperbolic systems.

In this work, we introduce a novel Hermite method to handle Maxwell's equations for nonlinear dispersive media. The proposed method achieves high-order accuracy and is free of any nonlinear algebraic solver, requiring solving instead small local linear systems for which the dimension is independent of the order. The implementation of order adaptive algorithms is straightforward in this setting, making the resulting p-adaptive Hermite method appealing for the simulations of soliton-like wave propagation.

In this work, we develop implicit rank-adaptive schemes for time-dependent matrix differential equations. The dynamic low rank approximation (DLRA) is a well-known technique to capture the dynamic low rank structure based on Dirac-Frenkel time-dependent variational principle. In recent years, it has attracted a lot of attention due to its wide applicability. Our schemes are inspired by the three-step procedure used in the rank adaptive version of the unconventional robust integrator (the so called BUG integrator) for DLRA. First, a prediction (basis update) step is made computing the approximate column and row spaces at the next time level. Second, a Galerkin evolution step is invoked using a base implicit solve for the small core matrix. Finally, a truncation is made according to a prescribed error threshold. Since the DLRA is evolving the differential equation projected on to the tangent space of the low rank manifold, the error estimate of the BUG integrator contains the tangent projection (modeling) error which cannot be easily controlled by mesh refinement. This can cause convergence issue for equations with cross terms. To address this issue, we propose a simple modification, consisting of merging the row and column spaces from the explicit step truncation method together with the BUG spaces in the prediction step. In addition, we propose an adaptive strategy where the BUG spaces are only computed if the residual for the solution obtained from the prediction space by explicit step truncation method, is too large. We prove stability and estimate the local truncation error of the schemes under assumptions. We benchmark the schemes in several tests, such as anisotropic diffusion, solid body rotation and the combination of the two, to show robust convergence properties.

Energy-conserving Hermite methods for solving Maxwell’s equations in dielectric and dispersive media are described and analyzed. In three space dimensions, methods of order 2m to 2m+2 require (m+1)^3 degrees-of-freedom per node for each field variable and can be explicitly marched in time with steps independent of m. We prove the stability for time steps limited only by domain-of-dependence requirements along with error estimates in a special semi-norm associated with the interpolation process. Numerical experiments are presented which demonstrate that Hermite methods of very high order enable the efficient simulation of the electromagnetic wave propagation over thousands of wavelengths.

We develop and analyze a new approach for simultaneously computing multiple solutions to the Helmholtz equation for different frequencies and different forcing functions. The new Multi-Frequency WaveHoltz (MFWH) algorithm is an extension of the original WaveHoltz method and both are based on time-filtering solutions to an associated wave equation. With MFWH, the different Helmholtz solutions are computed simultaneously by solving a single wave equation combined with multiple time filters. The MFWH algorithm defines a fixed-point iteration which can be accelerated with Krylov methods such as GMRES. The solution of the wave equation can be efficiently solved with either explicit time-stepping or implicit time-stepping using as few as five time-steps per period. When combined with an O(N) solver for the implicit equations, such a multigrid, the scheme has an O(N) solution cost when the frequencies are fixed and the number of grid points N increases. High-order accurate approximations in space are used together with second-order accurate approximations in time. We show how to remove time discretization errors so that the MFWH solutions converge to the corresponding solutions to the discretized Helmholtz problems. Numerical results are given using second-order accurate and fourth-accurate discretizations to confirm the convergence theory.

High order accurate Hermite methods for the wave equation on curvilinear domains are presented. Boundaries are treated using centered compatibility conditions rather than more standard one-sided approximations. Both first-order-in-time (FOT) and second-order-in-time (SOT) Hermite schemes are developed. Hermite methods use the solution and multiple derivatives as unknowns and achieve space-time orders of accuracy 2m-1 (FOT) and 2m (SOT) for methods using (m+1)^d degree of freedom per node in d dimensions. The compatibility boundary conditions (CBCs) are based on taking time derivatives of the boundary conditions and using the governing equations to replace the time derivatives with spatial derivatives. These resulting constraint equations augment the Hermite scheme on the boundary. The solvability of the equations resulting from the compatibility conditions are analyzed. Numerical examples demonstrate the accuracy and stability of the new schemes in two dimensions.