Time-dependent operators

This tutorial explains how to define time-dependent Hamiltonians – and more generally time-dependent operators – in Dynamiqs. There are currently four supported formats: constant operator, piecewise constant operator, constant operator modulated by a time-dependent factor, or arbitrary time-dependent operator defined by a function.

import dynamiqs as dq
import jax.numpy as jnp

The `TimeArray` type

In Dynamiqs, time-dependent operators are defined with TimeArray objects. These objects can be called at arbitrary times, and return the corresponding array at that time. For example to define the Hamiltonian $$ H_x(t)=\cos(2\pi t)\sigma_x $$

>>> f = lambda t: jnp.cos(2.0 * jnp.pi * t)
>>> Hx = dq.modulated(f, dq.sigmax())  # initialize a modulated time-array
>>> Hx(1.0)
Array([[0.+0.j, 1.+0.j],
       [1.+0.j, 0.+0.j]], dtype=complex64)
>>> Hx.shape
(2, 2)

Time-arrays support common arithmetic operations with scalars, regular arrays and other time-array objects. For example to define the Hamiltonian $$ H(t) = \sigma_z + 2 H_x(t) - \sin(\pi t) \sigma_y $$

>>> g = lambda t: jnp.sin(jnp.pi * t)
>>> Hy = dq.modulated(g, dq.sigmay())
>>> H = dq.sigmaz() + 2 * Hx - Hy
>>> H(1.0)
Array([[ 1.+0.j,  2.-0.j],
       [ 2.+0.j, -1.+0.j]], dtype=complex64)

Finally, time-arrays also support common utility functions, such as .conj(), or .reshape(). More details can be found in the TimeArray API page.

Defining a `TimeArray`

Constant operators

A constant operator is defined by $$ O(t) = O_0 $$ for any time $t$, where $O_0$ is an arbitrary operator. The most practical way to define constant operators is using array-like objects. They can also be instantiated as TimeArray instances using the dq.constant() function. For instance, to define the Pauli operator $H = \sigma_z$, you can use any of the following syntaxes:

DynamiqsNumPy arrayJAX arrayQuTiP QobjPython list

import dynamiqs as dq
H = dq.sigmaz()

import numpy as np
H = np.array([[1, 0], [0, -1]])

import jax.numpy as jnp
H = jnp.array([[1, 0], [0, -1]])

import qutip as qt
H = qt.sigmaz()

H = [[1, 0], [0, -1]]

Note

Common operators are available as utility functions, see the list of available operators in the Python API.

Piecewise constant operators

A piecewise constant (PWC) operator takes constant values over some time intervals. It is defined by $$ O(t) = \left(\sum_{k=0}^{N-1} c_k\; \Omega_{[t_k, t_{k+1}[}(t)\right) O_0 $$ where $c_k$ are constant values, $\Omega_{[t_k, t_{k+1}[}$ is the rectangular window function defined by $\Omega_{[t_a, t_b[}(t) = 1$ if $t \in [t_a, t_b[$ and $\Omega_{[t_a, t_b[}(t) = 0$ otherwise, and $O_0$ is a constant operator.

In Dynamiqs, PWC operators are defined by three array-like objects:

times: the time points $(t_0, \ldots, t_N)$ defining the boundaries of the time intervals, of shape (N+1,),
values: the constant values $(c_0, \ldots, c_{N-1})$ for each time interval, of shape (..., N),
array: the array defining the constant operator $O_0$, of shape (n, n).

To construct a PWC operator, these three arguments must be passed to the dq.pwc() function, which returns a TimeArray object. When called at some time $t$, this object then returns an array with shape (..., n, n). For example, let us define a PWC operator $H(t)$ with constant value $3\sigma_z$ for $t\in[0, 1[$ and $-2\sigma_z$ for $t\in[1, 2[$:

>>> times = [0.0, 1.0, 2.0]
>>> values = [3.0, -2.0]
>>> array = dq.sigmaz()
>>> H = dq.pwc(times, values, array)
>>> H
PWCTimeArray(shape=(2, 2), dtype=complex64)

The returned object can be called at different times:

$t = -1.0$$t=0.0$$t=0.5$$t=1.0$$t=1.5$$t=2.0$

>>> H(-1.0)
Array([[ 0.+0.j,  0.+0.j],
       [ 0.+0.j, -0.+0.j]], dtype=complex64)

>>> H(0.0)
Array([[ 3.+0.j,  0.+0.j],
       [ 0.+0.j, -3.+0.j]], dtype=complex64)

>>> H(0.5)
Array([[ 3.+0.j,  0.+0.j],
       [ 0.+0.j, -3.+0.j]], dtype=complex64)

>>> H(1.0)
Array([[-2.+0.j, -0.+0.j],
       [-0.+0.j,  2.-0.j]], dtype=complex64)

>>> H(1.5)
Array([[-2.+0.j, -0.+0.j],
       [-0.+0.j,  2.-0.j]], dtype=complex64)

>>> H(2.0)
Array([[ 0.+0.j,  0.+0.j],
       [ 0.+0.j, -0.+0.j]], dtype=complex64)

Note

The argument times must be sorted in ascending order, but does not need to be evenly spaced. When calling the resulting time-array object at time $t$, the returned array is the operator $c_k\ O_0$ corresponding to the interval $[t_k, t_{k+1}[$ in which the time $t$ falls. If $t$ does not belong to any time intervals, the returned array is null.

Batching PWC operators

The batching of the returned time-array is specified by values. For example, to define a PWC operator batched over a parameter $\theta$:

>>> thetas = jnp.linspace(0, 1.0, 11)  # (11,)
>>> times = [0.0, 1.0, 2.0]
>>> values = thetas[:, None] * jnp.array([3.0, -2.0])  # (11, 2)
>>> array = dq.sigmaz()
>>> H = dq.pwc(times, values, array)
>>> H.shape
(11, 2, 2)

Modulated operators

A modulated operator is defined by $$ O(t) = f(t) O_0 $$ where $f(t)$ is an time-dependent scalar. In Dynamiqs, modulated operators are defined by:

f: a Python function with signature f(t: float) -> Scalar | Array that returns the modulating factor $f(t)$ for any time $t$, as a scalar or an array of shape (...),
array: the array defining the constant operator $O_0$, of shape (n, n).

To construct a modulated operator, these two arguments must be passed to the dq.modulated() function, which returns a TimeArray object. When called at some time $t$, this object then returns an array with shape (..., n, n). For example, let us define the modulated operator $H(t)=\cos(2\pi t)\sigma_x$:

>>> f = lambda t: jnp.cos(2.0 * jnp.pi * t)
>>> H = dq.modulated(f, dq.sigmax())
>>> H
ModulatedTimeArray(shape=(2, 2), dtype=complex64)

The returned object can be called at different times:

$t = 0.5$$t=1.0$

>>> H(0.5)
Array([[-0.+0.j, -1.+0.j],
       [-1.+0.j, -0.+0.j]], dtype=complex64)

>>> H(1.0)
Array([[0.+0.j, 1.+0.j],
       [1.+0.j, 0.+0.j]], dtype=complex64)

Batching modulated operators

The batching of the returned time-array is specified by the array returned by f. For example, to define a modulated Hamiltonian $H(t)=\cos(\omega t)\sigma_x$ batched over the parameter $\omega$:

>>> omegas = jnp.linspace(0.0, 1.0, 11)  # (11,)
>>> f = lambda t: jnp.cos(omegas * t)
>>> H = dq.modulated(f, dq.sigmax())
>>> H.shape
(11, 2, 2)

Function with additional arguments

To define a modulated operator with a function that takes arguments other than time (extra *args and **kwargs), you can use functools.partial(). For example:

>>> import functools
>>> def pulse(t, omega, amplitude=1.0):
...     return amplitude * jnp.cos(omega * t)
>>> # create function with correct signature (t: float) -> Array
>>> f = functools.partial(pulse, omega=1.0, amplitude=5.0)
>>> H = dq.modulated(f, dq.sigmax())

Discontinuous function

If there is a discontinuous jump in the function values, you should use the optional argument discontinuity_ts to enforce adaptive step size solvers to stop at these times (i.e., right before, and right after the jump).

Arbitrary time-dependent operators

An arbitrary time-dependent operator is defined by $$ O(t) = f(t) $$ where $f(t)$ is a time-dependent operator. In Dynamiqs, arbitrary time-dependent operators are defined by:

f: a Python function with signature f(t: float) -> Array that returns the operator $f(t)$ for any time $t$, as an array of shape (..., n, n).

To construct an arbitrary time-dependent operator, pass this argument to the dq.timecallable() function, which returns a TimeArray object. This object then returns an array with shape (..., n, n) when called at any time $t$.

For example, let us define the arbitrary time-dependent operator $H(t)=\begin{pmatrix}t & 0\\0 & 1 - t\end{pmatrix}$:

>>> f = lambda t: jnp.array([[t, 0], [0, 1 - t]])
>>> H = dq.timecallable(f)
>>> H
CallableTimeArray(shape=(2, 2), dtype=float32)

The returned object can be called at different times:

$t = 0.5$$t=1.0$

>>> H(0.5)
Array([[0.5, 0. ],
       [0. , 0.5]], dtype=float32)

>>> H(1.0)
Array([[1., 0.],
       [0., 0.]], dtype=float32)

The function f must return a JAX array (not an array-like object!)

An error is raised if the function f does not return a JAX array. This error concerns any other array-like objects. This is enforced to avoid costly conversions at every time step of the numerical integration.

Batching arbitrary time-dependent operators

The batching of the returned time-array is specified by the array returned by f. For example, to define an arbitrary time-dependent operator batched over a parameter $\theta$:

>>> thetas = jnp.linspace(0, 1.0, 11)  # (11,)
>>> f = lambda t: thetas[:, None, None] * jnp.array([[t, 0], [0, 1 - t]])
>>> H = dq.timecallable(f)
>>> H.shape
(11, 2, 2)

Function with additional arguments

To define an arbitrary time-dependent operator with a function that takes arguments other than time (extra *args and **kwargs), you can use functools.partial(). For example:

>>> import functools
>>> def func(t, a, amplitude=1.0):
...     return amplitude * jnp.array([[t, a], [a, 1 - t]])
>>> # create function with correct signature (t: float) -> Array
>>> f = functools.partial(func, a=1.0, amplitude=5.0)
>>> H = dq.timecallable(f)

Discontinuous function

If there is a discontinuous jump in the function values, you should use the optional argument discontinuity_ts to enforce adaptive step size solvers to stop at these times (i.e., right before, and right after the jump).

Time-dependent operators

The TimeArray type

Defining a TimeArray

Constant operators

Piecewise constant operators

Modulated operators

Arbitrary time-dependent operators

The `TimeArray` type

Defining a `TimeArray`