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
inpropagate
/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 toQuantumPropagators
:ExpProp
– constructs $\op{U}$ explicitly and then applies it to $|Ψ⟩$Cheby
— expansion of $\op{U} |Ψ⟩$ into Chebychev polynomials, valid if $\op{H}$ has real eigenvaluesNewton
– 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 polynomials series expansions are more suitable for bigger systems and provide both efficiency and high precision (in general, the 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.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 settingpwc=true
, for comparison withmethod=Cheby
andmethod=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
— Methodusing 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 argumentH_dt
is obtained by constructing an operatorH
fromgenerator
via theevaluate
function and the multiplied with the time stepdt
for the current time interval. The propagation then simply multiplies the return value offunc
with the current stateconvert_state
: Type to which to temporarily convert states before multiplying the return value offunc
.convert_operator
: Type to which to convert the operatorH
before multiplying it withdt
and plugging the result intofunc
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.
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 (orH
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
— Methodusing 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 ingenerator
(seeget_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 thespecrange
functionspecrange_buffer
: An additional factor by which to enlarge the estimated spectral range returned byspecrange
, 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-zerocheck_normalization
: Check whether the Hamiltonian has been properly normalized, i.e., that the spectral range ofgenerator
has not been underestimated. This slowes down the propagation, but is advisable for novelgenerators
.uniform_dt_tolerance=1e-12
: How much the intervals oftlist
are allowed to vary while still being considered constant.specrange_kwargs
: All further keyword arguments are passed to thespecrange
function. Most notably, with the defaultspecrange_method=:auto
(orspecrange_method=:manual
), passingE_min
andE_max
allows to manually specify the spectral range ofgenerator
.
The time evolution operator of the piecewise-constant Schrödinger equation $|Ψ(t)⟩ = e^{-i Ĥ dt} |Ψ(0)⟩$ is evaluated by an expansion into Chebychev polynomials [1, 2]. This requires $Ĥ$ 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
Newton
The method should be loaded with
using QuantumPropagators: Newton
and then passed as method=Newton
to propagate
or init_prop
:
QuantumPropagators.init_prop
— Methodusing 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
The time evolution operator of the piecewise-constant Schrödinger equation $|Ψ(t)⟩ = e^{-i Ĥ dt} |Ψ(0)⟩$ is evaluated by an expansion into Newton polynomials [3, 5, 10]. Unlike for Chebychev polynomials, this expansion does not require $Ĥ$ to be Hermitian or to have a known spectral radius. This makes the Newton propagation applicable to open quantum systems, where $Ĥ$ 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
andmax_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
— Methodusing 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 defaultpwc=false
, for time-continuous controls.piecewise
: Currently equivalent topwc
, 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, theparameters
are determined automatically usingget_parameters
, respectivelydiscretize_on_midpoints
ifpwc=true
. If they are given manually, they must follow the exact same semantics. In particular, forpwc=false
, any parameters must alias the parameters in the controls, such that mutatingparameters
is automatically reflected inevaluate
. Theparameters
will be available as an attribute of theode_propagator
.alg
: The algorithm to use for the ODE Solver, see the list of solvers in the DifferentialEquations manual. The defaultTsit5()
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, sincealg
must be specified manually. Also, the "Output Control" is managed by theode_propagator
, so these options should not be used.
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