Time Propagation

The time-evolution of quantum systems is a fundamental problem in quantum mechanics that can be computationally very demanding, especially for large Hilbert and Liouville spaces. SimOS provides two basic routines to simulate the time-evolution of (open) quantum systems:

the evol() routine allows to simulate dynamics under time-independent Hamiltonians or Liouvillians
the prop() routine allows to simulate evolution of a quantum sysmtem under time-dependent Hamiltonian or Liouvillians. Both the coherent as well as the incoherent contributions may be time-dependent.

Time-independent evolution in SimOS

In the simplest case, the system Hamiltonian (or Liouvillian, if incoherent dynamics are present) is static, i.e. time-independent. In this case, the evol() routine can be utilized so simulate the time evolution of the quantum system. The method calculates the matrix exponential of the Hamiltonian/Liouvillian (i.e. the propagator \(U\)) and applies it to an initial state as explained in more detail in our QM-primer. If incoherent dynamics are included, a Lindblad form of the QME is utilized.

The function is called with the static system-Hamiltonian H , the time step (i.e. duration) of the evolution period t and, if incoherent dynamics are to be included in the simulation, a list of collapse operators c_ops (keyword-argument). If an initial state rho_0 is provided, the evolved state is returned. If the initial state is omitted, the function instead returns the propagator of the evolution. The latter is especially useful to calculate the evolution under the same static Hamiltonian for many time steps, since the costly matrix exponential only has to be evaluated once.

Note

To perform state rotations on a Bloch sphere, i.e. rotations with well-defined rotation angle and axis in the complex vector space of the quantum system, you may use the simos.propagation.rot routine instead. It functions identical to simos.propagation.evol, however it acccepts a rotation operator and rotation angle as input arguments and is therefore the more convenient interface for simple state rotations. In full analogy to the simos.propagation.evol method, state evolution is performed if an initial state is being provided while the propagator for the rotation is returned if the initial state is omitted. To visualize state rotations, you may also find our Bloch-sphere visualizer helpful.

Let us consider a simple example, a single electron spin 1/2 that is polarized along an external magnetic field \(B_0\). We want to perform a \(\pi/2\) pulse (i.e. a 90-degree rotation on the Bloch sphere) and evolve the state under a static Zeeman interaction with the external magnetic field for a single time step. For the 90-degree rotation, we utilize the rot() routine, for the evolution under the static Zeeman Hamiltonian we use evol().

# quantum system of a single electron spin
S  = {'name': 'S', 'val': 1/2}
s = sos.System([S])
# the initial state, polarized electron spin
rho0  = sos.state(s, 'S[0.5]')
# rotate state with a pi/2 pulse
rho0x = sos.rot(s.Sx, np.pi/2, rho0)
# evolve the state
dt = 1e-6 # time step [s]
B = 10e-3 # magnetic field [mT]
HZ = sos.ye*s.Sz # the Zeeman Hamiltonian
rhoxevol = sos.evol(HZ, dt, rho0x)

Time-dependent evolution in SimOS

Hamiltonian and collapse operators may become time-dependent under the application of time-varying control fields, e.g., magnetic or electric fields or shaped laser-excitation pulses. In most cases, the dynamic behavior can be efficiently parametrized with a limited set of control fields and collapse operators without any loss of generality. The time-varying Hamiltonian and collapse operators are formulated as

\[\hat{H}(t) = \hat{H}_0 + \sum_{i=1}^N c_i(t) \hat{H}_i\]

and

\[\hat{L}_k(t) = \hat{L}_{k,0} + \sum_{i=1}^N C_i(t) \hat{L}_{k,i}\]

with time-independent basis functions \(\hat{H}_i\) and \(\hat{L}_{k,i}\) and control amplitudes \(c_i(t)\) and \(C_i(t)\).

The prop() routine of SimOS accepts the time-independent parts \(\hat{H}_0\) (\(\hat{L}_{k,0}\)) as well as lists of time-varying control operators \(\hat{H}_i\) (\(\hat{L}_{k,i}\)) and their control amplitudes \(c_i(t)\) (\(C_i(t)\)) as input arguments. The time-dependence of the control amplitudes is assumed to be piecewise constant and control amplitudes are specified for series of discretized time-intervals \(dt\). A detailed explanation of all input parameters for prop() can be found in the syntax reference of the method.

Note

The prop() routine is only compatible with and thus implemented for our numerical backends and not for our symbolic Sympy backend. To propagate quantum systems in a symbolic matter, please utilize the evol() routine instead.

Engines

The development of computationally efficient integrator schemes for the numerical propagation of quantum systems is an active field of research. Besides classical ordinary differential equation (ODE) solvers, Euler-type integrators are common in the field of magnetic resonance simulations. Here, in full analogy to the static solutions presented in our QM-primer, the matrix exponential of the piecewise constant Hamiltonian or Liouvillian is evaluated separately for each time interval. The repeated computation of the matrix exponential is computationally expensive. However, a more efficient implementation can be achieved by using parallel-in-time integrators such as PARAMENT, introduced recently by some of use. We have further shown that Euler-type integrators can be easily converted into Magnus-type integrators that benefit from much quicker convergence.

The prop() routine serves as an interface to different so-called engines that perform the calculation and are readily selected via the engine keyword argument. You can thus select an integrator scheme that is optimal for your specific problem characteristics such as dimensionality and sparsity of the operators.

Engines for time propagation.
Engine	Compatible Backends	Description
`cpu`	All backends	CPU-based Euler-type or Magnus-type integrator (selected with keyword argument `Magnus`)
`parament`	All backends	GPU-based Magnus integrator. Utilizes the Parament package to perform parallelized time-propagation on a GPU, requires installed Parament and a GPU
`qutip`	Only available for qutip backend	Uses native QuTIP `mesolve` method, which is based on a Runge-Kutta solver
`RK45`	All backends.	Utilizes `scipy.integrate` 4th order Runge-Kutta

Note

All technical details of our Parament engine are available in the original publication. Esentially, Parament leverages the computation power of graphics processing units (GPUs) to perform the integration of all time steps in parallel.

Syntax Reference

evol(H, t, *rho, c_ops=[], proj=[], wallclock='global')[source]

Evolve state rho under static conditions for time t and return rho’ if rho is given.

Parameters:

H – Hamiltonian in angular frequencies for the evolution. Can be a quantum object or None if c_ops are given.
t – Evolution time.
*rho –
Input state (optional)

Returns:

Evolved state. If rho is not specified, the evolution operator is returned. In cases where c_ops are given, the vector superoperator is returned.

Keyword Arguments:

c_ops (list) –
List of collapse operators. If c_ops is given, the evolution is following the Lindblad master equation.
wallclock –
If ‘global’, the global wallclock is used to track the time. If None, no time tracking is performed. If a wallclock object is given, it is used to track the time.
proj (list) –
List of projectors for space reduction. Currently not implemented, reserved for future use.

prop(H0, dt, *rho, c_ops=[], H1=None, carr1=None, c_ops2=[], carr2=None, proj=[], engine='cpu', wallclock='global', **kwargs)[source]

Evolve state for time-dependent conditions for given parametrization and return rho’ if rho is given.

Parameters:

H0 – Static background Hamiltonian in angular frequencies for the evolution. Can be a quantum object or None. Note, at least one of H0, H1, c_ops or c_ops2 must be given.
dt – Time step for the modulation array(s) for the evolution.
*rho –
Input state (optional)

Returns:

Evolved state. If rho is not specified, the evolution operator is returned. In cases where c_ops are given, the vector superoperator is returned.

Keyword Arguments:

c_ops (list) –
List of collapse operators. If c_ops is given, the evolution is following the Lindblad master equation.
H1 (list) –
List of time-dependent Hamiltonians. If H1 is given, carr1 must be given as well. The number of Hamiltonians in H1 and the number of modulation arrays in carr1 must be equal.
carr1 (list) –
List of coefficients for the time-dependent Hamiltonians. If carr1 is given, H1 must be given as well. The number of arrays in carr1 and the number of Hamiltonians in H1 must be equal. Discretization of the time-dependent Hamiltonians is done by the time step dt.
c_ops2 (list) –
List of time-dependent collapse operators. If c_ops2 is given, carr2 must be given as well. The number of collapse operators in c_ops2 and the number of modulation arrays in carr2 must be equal.
carr2 (list) –
List of coefficients for the time-dependent collapse operators. If carr2 is given, c_ops2 must be given as well. The number of arrays in carr2 and the number of collapse operators in c_ops2 must be equal. Discretization of the time-dependent collapse operators is done by the time step dt.
engine (str) –
Engine to be used for the evolution. Can be ‘cpu’, ‘qutip’ or ‘parament’. If ‘qutip’ is chosen, mesolve from qutip is used for the evolution. This only works if the qutip backend is used. If ‘parament’ is chosen, the parament library is used for the evolution. This only works if the parament backend is installed and a GPU is available to use.
wallclock –
If ‘global’, the global wallclock is used to track the time. If None, no time tracking is performed. If a wallclock object is given, it is used to track the time.
proj (list) –
List of projectors for space reduction. Currently not implemented, reserved for future use.
kwargs –
Additional keyword arguments for the evolution. Will be passed to the used engine. For the qutip engine, e.g., this can be the options for the mesolve function.

This operation corresponds to an infinitely fast pulse of the spin operator along the axis defined by the operator op by an angle phi.