Graph.density_matrix_evolution_pwc(initial_density_matrix, hamiltonian, lindblad_terms, sample_times=None, error_tolerance=1e-06, *, name=None)

Calculate the state evolution of an open system described by the GKS–Lindblad master equation.

The controls that you provide to this function have to be in piecewise-constant format. If your controls are smooth sampleable tensor-valued functions (STFs), you have to discretize them with discretize_stf before passing them to this function. You may need to increase the number of segments that you choose for the discretization depending on the sizes of oscillations in the smooth controls.

By default, this function computes an approximate piecewise-constant solution for the consideration of efficiency, with the accuracy controlled by the parameter error_tolerance. If your system is small, you can set the error_tolerance to None to obtain an exact solution. However, note that the exact method could be slow and memory intensive when applied to large systems.

  • initial_density_matrix (np.ndarray or Tensor) – A 2D array of the shape (D, D) representing the initial density matrix of the system, \(\rho_{\rm s}\). You can also pass a batch of density matrices and the input data shape must be (B, D, D) where B is the batch dimension.

  • hamiltonian (Pwc or SparsePwc) – A piecewise-constant function representing the effective system Hamiltonian, \(H_{\rm s}(t)\), for the entire evolution duration. If you pass any Lindblad operator as a dense array, the Hamiltonian will get converted to a (dense) Pwc.

  • lindblad_terms (list[tuple[float, np.ndarray or Tensor or scipy.sparse.spmatrix]]) – A list of pairs, \((\gamma_j, L_j)\), representing the positive decay rate \(\gamma_j\) and the Lindblad operator \(L_j\) for each coupling channel \(j\). You must provide at least one Lindblad term. If you pass the Hamiltonian as a Pwc, the operators will get converted to dense operators.

  • sample_times (np.ndarray, optional) – A 1D array of length \(T\) specifying the times \(\{t_i\}\) at which this function calculates system states. Must be ordered and contain at least one element. Note that increasing the density of sample times does not affect the computation precision of this function.

  • error_tolerance (float, optional) – Defaults to 1e-6. This option enables an approximate method to solve the master equation, meaning the 2-norm of the difference between the propagated state and the exact solution at the final time (and at each sample time if passed) is within the error tolerance. Note that, if set, this value must be smaller than 1e-2 (inclusive). However, setting it to a too small value (for example below 1e-12) might result in slower computation, but would not further improve the precision, since the dominating error in that case is due to floating point error. You can also explicitly set this option to None to find the exact piecewise-constant solution. Note that using the exact solution can be very computationally demanding in calculations involving a large Hilbert space or a large number of segments.

  • name (str, optional) – The name of the node.


System states at sample times. The shape of the return value is (D, D) or (T, D, D), depending on whether you provide sample times. Otherwise, the shape is (B, T, D, D) if you provide a batch of initial states.

Return type


See also


Discretize an Stf into a Pwc.


Create SparsePwc operators.


Corresponding operation for coherent evolution.


Under Markovian approximation, the dynamics of an open quantum system can be described by the GKS–Lindblad master equation 1 2

\[\frac{{\rm d}\rho_{\rm s}(t)}{{\rm d}t} = -i [H_{\rm s}(t), \rho_{\rm s}(t)] + \sum_j \gamma_j {\mathcal D}[L_j] \rho_{\rm s}(t) ,\]

where \({\mathcal D}\) is a superoperator describing the decoherent process in the system evolution and defined as

\[{\mathcal D}[X]\rho := X \rho X^\dagger - \frac{1}{2}\left( X^\dagger X \rho + \rho X^\dagger X \right)\]

for any system operator \(X\).

This function uses sparse matrix multiplication when the Hamiltonian is passed as a SparsePwc and the Lindblad operators as sparse matrices. This leads to more efficient calculations when they involve large operators that are relatively sparse (contain mostly zeros). In this case, the initial density matrix is still a densely represented array or tensor.



V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).


G. Lindblad, Commun. Math. Phys. 48, 119 (1976).


Simulate a trivial decay process for a single qubit described by the following master equation \(\dot{\rho} = -i[\sigma_z / 2, \, \rho] + \mathcal{D}[\sigma_-]\rho\).

>>> duration = 20
>>> initial_density_matrix = np.array([[0, 0], [0, 1]])
>>> hamiltonian = graph.constant_pwc_operator(
...     duration=duration, operator=graph.pauli_matrix("Z") / 2
... )
>>> lindblad_terms = [(1, graph.pauli_matrix("M"))]
>>> graph.density_matrix_evolution_pwc(
...     initial_density_matrix, hamiltonian, lindblad_terms, name="decay"
... )
<Tensor: name="decay", operation_name="density_matrix_evolution_pwc", shape=(2, 2)>
>>> result = qctrl.functions.calculate_graph(graph=graph, output_node_names=["decay"])
>>> result.output["decay"]["value"]
array([[9.99999998e-01+0.j, 0.00000000e+00+0.j],
       [0.00000000e+00+0.j, 2.06115362e-09+0.j]])

See more examples in the How to simulate open system dynamics and How to simulate large open system dynamics user guides.