Graph.steady_state(hamiltonian, lindblad_terms, method='QR', *, name=None)

Find the steady state of a time-independent open quantum system.

The Hamiltonian and Lindblad operators that you provide have to give rise to a unique steady state.

  • hamiltonian (Tensor or spmatrix) – A 2D array of shape (D, D) representing the time-independent Hamiltonian of the system, \(H_{\rm s}\).

  • 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.

  • method (str, optional) – The method used to find the steady state. Must be one of “QR”, “EIGEN_SPARSE”, or “EIGEN_DENSE”. The “QR” method obtains the steady state through a QR decomposition and is suitable for small quantum systems with dense representation. The “EIGEN_SPARSE” and “EIGEN_DENSE” methods find the steady state as the eigenvector whose eigenvalue is closest to zero, using either a sparse or a dense representation of the generator. Defaults to “QR”.

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


The density matrix representing the steady state of the system.

Return type:



This function currently does not support calculating the gradient with respect to its inputs. Therefore, it cannot be used in a graph for a calculate_optimization or calculate_stochastic_optimization call, which will raise a RuntimeError. Please use gradient-free optimization if you want to perform an optimization task with this function. You can learn more about it in the How to optimize controls using gradient-free optimization user guide.

See also


State evolution of an 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],

\[\frac{{\rm d}\rho_{\rm s}(t)}{{\rm d}t} = {\mathcal L} (\rho_{\rm s}(t)) ,\]

where the Lindblad superoperator \({\mathcal L}\) is defined as

\[{\mathcal L} (\rho_{\rm s}(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 computes the steady state of \({\mathcal L}\) by solving

\[\frac{{\rm d}\rho_{\rm s}(t)}{{\rm d}t} = 0 .\]

The function assumes that \(H_{\rm s}\) is time independent and that the dynamics generated by \({\mathcal L}\) give rise to a unique steady state. That is, the generated quantum dynamical map has to be ergodic [3].



Compute the steady state of the single qubit open system dynamics according to the Hamiltonian \(H=\omega\sigma_z\) and the single Lindblad operator \(L=\sigma_-\).

>>> omega = 0.8
>>> gamma = 0.5
>>> hamiltonian = omega * graph.pauli_matrix("Z")
>>> lindblad_terms = [(gamma, graph.pauli_matrix("M"))]
>>> graph.steady_state(hamiltonian, lindblad_terms, name="steady_state")
<Tensor: name="steady_state", operation_name="steady_state", shape=(2, 2)>
>>> result = qctrl.functions.calculate_graph(graph=graph, output_node_names=["steady_state"])
>>> result.output["steady_state"]["value"]
array([[1.+0.j 0.-0.j]
       [0.-0.j 0.-0.j]])