Skip to content

dq.steadystate 

steadystate(
    H: QArrayLike | ConstantTimeQArray,
    jump_ops: list[QArrayLike | ConstantTimeQArray],
    *,
    exp_ops: list[QArrayLike] | None = None,
    method: Method = TimeEvolution(),
    gradient: Gradient | None = None,
    cartesian_batching: bool = True,
    save_extra: Callable[[QArray], PyTree] | None = None,
    vectorized: bool = False,
    assume_hermitian: bool = True
) -> SteadyStateResult

Compute the steady state of the Lindblad master equation.

This function computes a steady state \(\rho_\infty\) of the Lindblad master equation (with \(\hbar=1\)), defined by \(\mathcal{L}(\rho_\infty) = 0\) with $$ \mathcal{L}(\rho) = -i[H, \rho] + \sum_{k=1}^N \left( L_k \rho L_k^\dag - \frac{1}{2} L_k^\dag L_k \rho - \frac{1}{2} \rho L_k^\dag L_k \right), $$ where \(H\) is the system's Hamiltonian and \(\{L_k\}\) is a collection of jump operators, both time-independent.

Parameters:

  • H (qarray-like or constant timeqarray of shape (...H, n, n)) –

    Hamiltonian.

  • exp_ops (list of qarray-like, each of shape (n, n)) –

    List of operators for which the expectation value is computed.

  • method –

    Method for the computation. Defaults to dq.method.TimeEvolution (supported: TimeEvolution).

  • gradient –

    Algorithm used to compute the gradient. The default is method-dependent, refer to the documentation of the chosen method for more details.

Returns:

  • dq.SteadyStateResult object holding the result of the steady state computation. Use result.state to access the steady state, result.expects to access its expectation values, and result.infos.converged to check that the method converged.

    Detailed result API 
    dq.SteadyStateResult

    Attributes:

    • state (qarray of shape (..., n, n)) - Steady state, or last state reached if the method did not converge.
    • expects (array of shape (..., len(exp_ops)) or None) - Expectation values in the steady state, if specified by exp_ops.
    • extra (PyTree or None) - Extra data saved with save_extra() if specified.
    • infos (PyTree) - Information on the resolution: infos.converged (bool array of shape (...)) is True if the convergence criterion of the method was met, infos.time (array of shape (...)) is the time at which the method stopped, and infos.steps holds the number of steps.
    • method (Method) - Method used.
    • gradient (Gradient) - Gradient used.
    • options (Options) - Options used.

Other Parameters:

  • cartesian_batching –

    If True, batched arguments are treated as separated batch dimensions, otherwise the batching is performed over a single shared batch dimension. Defaults to True.

  • save_extra –

    Function with signature f(QArray) -> PyTree that takes the steady state as input and returns a PyTree, accessible in result.extra. Defaults to None.

  • vectorized –

    If True, the master equation is solved by vectorizing the density matrix and Liouvillian, see dq.mesolve(). Defaults to False.

  • assume_hermitian –

    If True, the initial density matrix of the method is assumed to be Hermitian, see dq.mesolve(). Defaults to True.

Examples:

import dynamiqs as dq

n = 16
a = dq.destroy(n)

H = a.dag() @ a + 0.5 * (a + a.dag())
jump_ops = [a]

result = dq.steadystate(H, jump_ops, exp_ops=[a])
print(result)
Output
==== SteadyStateResult ====
Method  : TimeEvolution
Infos   : converged at t=18.7 | 180 steps (164 accepted, 16 rejected)
State   : QArray complex64 (16, 16) | 2.0 Kb
Expects : Array complex64 (1,) | 0.0 Kb

Advanced use-cases

Running multiple computations concurrently

The Hamiltonian H and the jump operators jump_ops can be batched to compute multiple steady states concurrently. All other arguments are common to every batch. The resulting states and expectation values are batched according to the leading dimensions of H and jump_ops. The behaviour depends on the value of the cartesian_batching option.

The results leading dimensions are 

... = ...H, ...L0, ...L1, (...)

For example if:

  • H has shape (2, 3, n, n),
  • jump_ops = [L0, L1] has shape [(4, 5, n, n), (6, n, n)],

then result.state has shape (2, 3, 4, 5, 6, n, n).

The results leading dimensions are 

... = ...H = ...L0 = ...L1 = (...)  # (once broadcast)

For example if:

  • H has shape (2, 3, n, n),
  • jump_ops = [L0, L1] has shape [(3, n, n), (2, 1, n, n)],

then result.state has shape (2, 3, n, n).

See the Batching simulations tutorial for more details.