Propagation Methods

\[\gdef\op#1{\hat{#1}} \gdef\ket#1{\vert{#1}\rangle} \gdef\Liouvillian{\mathcal{L}} \gdef\Re{\operatorname{Re}} \gdef\Im{\operatorname{Im}}\]

As discussed in the Overview, time propagation can be implemented in one of two ways:

By analytically solving the equation of motion and numerically evaluating the application time evolution operator.
We consider this especially in the piecewise-constant case (pwc=true in propagate/init_prop), which is required for the traditional optimization methods GRAPE and Krotov. In these propagations, the time-dependent generator $\op{H}(t)$ is evaluated to a constant operator $\op{H}$ on each interval of the time grid. The analytical solution to the Schrödinger or Liouville equation is well known, and propagation step simply has to evaluate the application of the time evolution operator $\op{U} = \exp[-i \op{H} dt]$ to the state $|Ψ⟩$. The following methods are built in to QuantumPropagators:
- ExpProp – constructs $\op{U}$ explicitly and then applies it to $|Ψ⟩$
- Cheby — expansion of $\op{U} |Ψ⟩$ into Chebychev polynomials, valid if $\op{H}$ has real eigenvalues
- Newton – expansion of $\op{U} |Ψ⟩$ into Newton polynomials, valid if $\op{H}$ has complex eigenvalues (non-Hermitian Hamiltonian, Liouvillian)
The ExpProp method is generally not numerically efficient, but works well for small system for for debugging. The two "core" methods based on a polynomial series expansions are more suitable for bigger systems and provide both efficiency and high precision (in general, the series is truncated as soon as some desired precision is reached, which is machine precision by default).
Note that this high precision is within the piecewise-constant approximation. The discretization itself may introduce a non-negligible error compared to the time-continuous dynamics. There is tradeoff: A smaller dt decreases the discretization error, but the polynomial expansions are more effective with larger time steps.
There is also support for external packages to efficiently evaluate the application of the piecewise-constant time evolution operator:
- ExponentialUtilities – applies $\op{U} |Ψ⟩$ using Krylov methods without explicitly forming $\op{U}$
By solving the equation of motion explicitly with an ODE solver.
We support the use of any of the ODE solvers in OrdinaryDiffEq.jl:
- OrdinaryDiffEq – solve the equation of motion as an ODE
The main benefit of using an ODE solver is that the generator $\op{H}(t)$ can be treated as time-continuous, and thus avoid the time discretization error. While this is not compatible with traditional optimal control method like GRAPE and Krotov, it is suitable for control methods for tuning analytical control parameters [6–8].
The method=OrdinaryDiffEq is also available in a piecewise-constant mode by setting pwc=true, for comparison with method=Cheby and method=Newton.

`ExpProp`

The method should be loaded with

using QuantumPropagators: ExpProp

and the passed as method=ExpProp to propagate or init_prop:

QuantumPropagators.init_prop — Method

using QuantumPropagators: ExpProp

exp_propagator = init_prop(
    state,
    generator,
    tlist;
    method=ExpProp,
    inplace=QuantumPropagators.Interfaces.supports_inplace(state),
    backward=false,
    verbose=false,
    parameters=nothing,
    func=(H_dt -> exp(-1im * H_dt))
    convert_state=_exp_prop_convert_state(state),
    convert_operator=_exp_prop_convert_operator(generator),
    _...
)

initializes an ExpPropagator.

Method-specific keyword arguments

func: The function to evaluate. The argument H_dt is obtained by constructing an operator H from generator via the evaluate function and the multiplied with the time step dt for the current time interval. The propagation then simply multiplies the return value of func with the current state
convert_state: Type to which to temporarily convert states before multiplying the return value of func.
convert_operator: Type to which to convert the operator H before multiplying it with dt and plugging the result into func

The convert_state and convert_operator parameters are useful for when the generator and or state are unusual data structures for which the relevant methods to calculate func are not defined. Often, it is easier to temporarily convert them to standard complex matrices and vectors than to implement the missing methods.

source

Advantages

Simple: no knobs to turn
"Exact" to machine precision (within the piecewise constant approximation)
Does not require any special properties or knowledge of the dynamical generator
Efficient for small systems

Disadvantages

Bad numerical scaling with the Hilbert space dimension
Method for exp(-1im * H * dt) must be defined (or H must be convertible to a type that can be exponentiated)

When to use

Small Hilbert space dimension (<10)
Comparing against another propagator

`Cheby`

The method should be loaded with

using QuantumPropagators: Cheby

and then passed as method=Cheby to propagate or init_prop:

QuantumPropagators.init_prop — Method

using QuantumPropagators: Cheby

cheby_propagator = init_prop(
    state,
    generator,
    tlist;
    method=Cheby,
    inplace=QuantumPropagators.Interfaces.supports_inplace(state),
    backward=false,
    verbose=false,
    parameters=nothing,
    control_ranges=nothing,
    specrange_method=:auto,
    specrange_buffer=0.01,
    cheby_coeffs_limit=1e-12,
    check_normalization=false,
    specrange_kwargs...
)

initializes a ChebyPropagator.

Method-specific keyword arguments

control_ranges: a dict the maps the controls in generator (see get_controls) to a tuple of min/max values. The Chebychev coefficients will be calculated based on a spectral envelope that assumes that each control can take arbitrary values within the min/max range. If not given, the ranges are determined automatically. Specifying manual control ranges can be useful when the the control amplitudes (parameters) may change during the propagation, e.g. in a sequential-update control scheme.
specrange_method: Method to pass to the specrange function
specrange_buffer: An additional factor by which to enlarge the estimated spectral range returned by specrange, in order to ensure that Chebychev coefficients are based on an overestimation of the spectral range.
cheby_coeffs_limit: The maximum magnitude of Chebychev coefficients that should be treated as non-zero
check_normalization: Check whether the Hamiltonian has been properly normalized, i.e., that the spectral range of generator has not been underestimated. This slowes down the propagation, but is advisable for novel generators.
uniform_dt_tolerance=1e-12: How much the intervals of tlist are allowed to vary while still being considered constant.
specrange_kwargs: All further keyword arguments are passed to the specrange function. Most notably, with the default specrange_method=:auto (or specrange_method=:manual), passing E_min and E_max allows to manually specify the spectral range of generator.

source

The time evolution operator of the piecewise-constant Schrödinger equation $|Ψ(t)⟩ = e^{-i Ĥ dt} |Ψ(0)⟩$ is evaluated by an expansion into Chebychev polynomials [1, 2]. This requires $Ĥ$ to be Hermitian (have real eigenvalues) and to have a known spectral range, so that it can be normalized to the domain $[-1, 1]$ on which the Chebychev polynomials are defined.

See [9, Chapter 3.2.1] for a detailed description of the method.

Advantages

Very efficient for high precision and moderately large time steps

Disadvantages

Only valid for Hermitian operators
Must be able to estimate the spectral envelope

When to use

Closed quantum systems with piecewise constant dynamics

`ExponentialUtilities`

The method requires that the ExponentialUtilities.jl package is loaded

using ExponentialUtilities

and then passed as method=ExponentialUtilities to propagate or init_prop:

QuantumPropagators.init_prop — Method

using ExponentialUtilities

expv_propagator = init_prop(
    state,
    generator,
    tlist;
    method=ExponentialUtilities,
    inplace=QuantumPropagators.Interfaces.supports_inplace(state),
    backward=false,
    verbose=false,
    parameters=nothing,
    expv_kwargs=(;),
    _...
)

initializes an ExponentialUtilitiesPropagator.

Method-specific keyword arguments

expv_kwargs: NamedTuple of keyword arguments forwarded to the underlying ExponentialUtilities routines. For in-place propagation, these are passed to both arnoldi! (Krylov subspace construction) and expv! (exponentiation). For out-of-place propagation, they are passed to expv, which forwards them internally to arnoldi. The most useful keyword arguments are:
- m::Int: Maximum Krylov subspace dimension. Defaults to min(30, size(A, 1)). Larger values improve accuracy but increase memory and computation cost per step.
- ishermitian::Bool: Whether the operator is Hermitian. Defaults to LinearAlgebra.ishermitian(A). When true, the cheaper Lanczos iteration is used instead of the full Arnoldi process. Note that the propagator passes -𝕚 dt as the time argument, so a Hermitian generator $Ĥ$ still results in a non-Hermitian product $-𝕚 Ĥ dt$; however, ExponentialUtilities handles the complex time scaling internally.
- tol::Real: Tolerance for happy breakdown. Defaults to 1e-7. The Arnoldi iteration stops early when the norm of the next Krylov vector drops below tol * opnorm.
- opnorm: Operator norm of A. By default, this is computed automatically. Supplying it avoids redundant norm computations.
- iop::Int: Incomplete orthogonalization procedure length. Defaults to 0 (full orthogonalization). A positive value limits the number of previous Krylov vectors used for orthogonalization, reducing cost at the expense of numerical stability.
- mode::Symbol: Termination strategy for expv. Either :happy_breakdown (default) or :error_estimate. In :happy_breakdown mode, the iteration relies on early termination when the Krylov subspace captures the relevant dynamics. In :error_estimate mode, an adaptive step-size strategy is used with additional tolerance parameters rtol (relative tolerance, defaults to √tol).

source

This method evaluates $\exp(-i \op{H} dt) |Ψ⟩$ via a Krylov expv algorithm without explicitly forming the matrix exponential. It builds a Krylov subspace (via Arnoldi or Lanczos iteration) and then computes the action of the matrix exponential on the state within that subspace. This is often a good fit for larger systems or matrix-free operators where direct matrix exponentiation is too costly.

The propagator requires that states and operators satisfy the AbstractArray interface (specifically, supports_vector_interface for states and supports_matrix_interface for operators).

Advantages

Avoids explicit construction of $\op{U}$
Works with matrix-free operators that support mul!
Good scaling for large sparse systems
Supports both Hermitian (Lanczos) and non-Hermitian (Arnoldi) generators

Disadvantages

Requires ExponentialUtilities.jl
Performance depends on Krylov subspace dimension and operator structure
Requires operators and states to implement an array interface, see QuantumPropagators.Interfaces.supports_matrix_interface and QuantumPropagators.Interfaces.supports_vector_interface, respectively

When to use

Large, sparse, or matrix-free generators
Systems where $\op{U}$ is too expensive to form explicitly

`Newton`

The method should be loaded with

using QuantumPropagators: Newton

and then passed as method=Newton to propagate or init_prop:

QuantumPropagators.init_prop — Method

using QuantumPropagators: Newton

newton_propagator = init_prop(
    state,
    generator,
    tlist;
    method=Newton,
    inplace=QuantumPropagators.Interfaces.supports_inplace(state),
    backward=false,
    verbose=false,
    parameters=nothing,
    m_max=10,
    func=(z -> exp(-1im * z)),
    norm_min=1e-14,
    relerr=1e-12,
    max_restarts=50,
    _...
)

initializes a NewtonPropagator.

Method-specific keyword arguments

m_max: maximum Krylov dimension, cf. NewtonWrk
func, norm_min, relerr, max_restarts: parameter to pass to newton!

source

The time evolution operator of the piecewise-constant Schrödinger equation $|Ψ(t)⟩ = e^{-i Ĥ dt} |Ψ(0)⟩$ is evaluated by an expansion into Newton polynomials [3, 5, 10]. Unlike for Chebychev polynomials, this expansion does not require $Ĥ$ to be Hermitian or to have a known spectral radius. This makes the Newton propagation applicable to open quantum systems, where $Ĥ$ is replaced by a Liouvillian to calculate the time evolution of the density matrix.

See [9, Chapter 3.2.2] for a detailed description of the method.

Advantages

Reasonably efficient for high precision and moderately large time steps
Spectral radius does not need to be known

Disadvantages

Need to choose m_max and max_restarts well for good performance.

When to use

Open quantum systems with piecewise constant dynamics

`OrdinaryDiffEq`

The method requires that the OrdinaryDiffEq package is loaded

using OrdinaryDiffEq

Equivalently, the more general DifferentialEquations.jl package can be used.

using DifferentialEquations

There is no difference between these two options: OrdinaryDiffEq is just a smaller dependency, but DifferentialEquations may be preferred if the large DifferentialEquations framework is required for the project.

In any case, the loaded package to propagate or init_prop via the method keyword argument:

QuantumPropagators.init_prop — Method

using OrdinaryDiffEq  # or: `using DifferentialEquations`

ode_propagator = init_prop(
    state,
    generator,
    tlist;
    method=OrdinaryDiffEq,  # or: `method=DifferentialEquations`
    inplace=QuantumPropagators.Interfaces.supports_inplace(state),
    backward=false,
    verbose=false,
    parameters=nothing,
    piecewise=false,
    pwc=false,
    alg=OrdinaryDiffEq.Tsit5(),
    solver_options...
)

initializes a propagator that uses an ODE solver from the OrdinaryDiffEq.jl package as a backend.

By default, the resulting propagator is for time-continuous controls that can be evaluated with evaluate(control, t) for any t in the range of tlist[begin] to tlist[end]. The controls may be parametrized, see get_parameters. Any parameters will be available in the parameters attribute of the resulting ode_propagator, as a dictionary mapping controls to a vector of parameter values. Mutating ode_propagator.parameters[control] will be reflected in any subsequent call to prop_step!.

If pwc=true (or, equivalently piecewise=true), all controls will be discretized with discretize_on_midpoints and the propagation will be for piecewise constant dynamics. The resulting ode_propagator will be an instance of PWCPropagator, with the corresponding semantics. In particular, the ode_propatator.parameters will be a mapping of controls to discretized pulse values, not the analytical parameters obtained with get_parameters(control) as in the default case.

Internally, the generator will be wrapped with QuantumPropagators.ode_function. The resulting function f will be called internally as f(ϕ, Ψ, vals_dict, t) or f(Ψ, vals_dict, t) depending on the inplace keyword argument.

Method-specific keyword arguments

pwc: Whether to propagate for piecewise-constant controls or, with the default pwc=false, for time-continuous controls.
piecewise: Currently equivalent to pwc, but future version may change this to allow for other piecewise (e.g., piecewise-linear) controls.
parameters: If given, a mapping of controls to parameter values (pwc=false) or pulse values on the intervals of the time grid (pwc=true). By default, the parameters are determined automatically using get_parameters, respectively discretize_on_midpoints if pwc=true. If they are given manually, they must follow the exact same semantics. In particular, for pwc=false, any parameters must alias the parameters in the controls, such that mutating parameters is automatically reflected in evaluate. The parameters will be available as an attribute of the ode_propagator.
alg: The algorithm to use for the ODE Solver, see the list of solvers in the DifferentialEquations manual. The default Tsit5() method is the recommended choice for non-stiff problems.
solver_options: All other keyword arguments are passed to the ODE solver, see the list of Solve Keyword Arguments in the DifferentialEquations manual. Note that the options for "Default Algorithm Hinting" do not apply, since alg must be specified manually. Also, the "Output Control" is managed by the ode_propagator, so these options should not be used.

source

Advantages

Suitable for time-continuous dynamics
The full power of the DifferentialEquations.jl ecosystem
Efficient for moderate precisions

Disadvantages

Less efficient for piecewise-constant dynamics, and thus less suitable of PWC control methods

When to use

Time-continuous dynamics