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, HsH_{\rm s}
  • 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
  • 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 run_optimization or run_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.


Graph.density_matrix_evolution_pwc : State evolution of an open quantum system.

Graph.jump_trajectory_evolution_pwc : Trajectory-based 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,

dρs(t)dt=L(ρs(t)), \frac{{\rm d}\rho_{\rm s}(t)}{{\rm d}t} = {\mathcal L} (\rho_{\rm s}(t)) ,

where the Lindblad superoperator L{\mathcal L}

L(ρs(t))=i[Hs(t),ρs(t)]+jγjD[Lj]ρs(t), {\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 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 computes the steady state of L{\mathcal L}

dρs(t)dt=0. \frac{{\rm d}\rho_{\rm s}(t)}{{\rm d}t} = 0 .

The function assumes that HsH_{\rm s} is time independent and that the dynamics generated by L{\mathcal L}3.


[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).

[3] D. Burgarth, G. Chiribella, V. Giovannetti, P. Perinotti, and K. Yuasa, New J. Phys. 15 073045 (2013).


Compute the steady state of the single qubit open system dynamics according to the Hamiltonian H=ωσzH=\omega\sigma_z and the single Lindblad operator L=σ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 = bo.execute_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]])

Was this useful?