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.

Note that when using the exact method, both hamiltonian and lindblad_terms are converted to the dense representation, regardless of their original formats. This means that the computation can be slow and memory intensive when applied to large systems.

When using the approximate method, the sparse representation is used internally if hamiltonian is a SparsePwc, otherwise the dense representation is used.


  • initial_density_matrix (np.ndarray or Tensor) – A 2D array of the shape (D, D) representing the initial density matrix of the system, ρs\rho_{\rm s}(B, D, D) where B is the batch dimension.
  • hamiltonian (Pwc or SparsePwc) – A piecewise-constant function representing the effective system Hamiltonian, Hs(t)H_{\rm s}(t)
  • lindblad_terms (list [ tuple [ float , np.ndarray or Tensor or scipy.sparse.spmatrix ] ]) – A list of pairs, (γj,Lj)(\gamma_j, L_j), representing the positive decay rate γj\gamma_j and the Lindblad operator LjL_j for each coupling channel jj
  • sample_times (np.ndarray or None , optional) – A 1D array like object of length TT specifying the times {ti}\{t_i\}
  • error_tolerance (float or None , 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 or None , optional) – The name of the node.


The time-evolved density matrix, with shape (D, D) or (T, D, D), depending on whether you provided sample times. If you provide a batch of initial states, the shape is (B, T, D, D) or (B, D, D).

Return type



Graph.jump_trajectory_evolution_pwc : Trajectory-based state evolution of an open quantum system.

Graph.steady_state : Compute the steady state of open quantum system.


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

dρs(t)dt=i[Hs(t),ρs(t)]+jγjD[Lj]ρs(t), \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 D{\mathcal D}

D[X]ρ:=XρX12(XXρ+ρXX) {\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 XX

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.


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

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


Simulate a trivial decay process for a single qubit described by the following master equation ρ˙=i[σz/2,,ρ]+D[σ]ρ\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 = bo.execute_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.

Was this useful?