QuantumOptics.jl implements two different ways to calculate steady states. The first one is to perform a time evolution according to a master equation until an adequate accuracy is reached

Most of the arguments are identical to the ones used for a standard Master time evolution. There are a few differences, however. First of all it is no longer required to give an initial state, but rather optional. If no initial state is given it defaults to a density matrix with a single 1 in the first entry. Note, that while for a cavity this corresponds to a Fock state with zero photons, it might be undesirable to use the default value when working with atoms (spin particles), since this state corresponds to the (highest) excited state.

Additionally, one can explicitly specify the trace distance criterion tol, which is checked in every time step. The time evolution stops, i.e. steady-state is reached, when the trace distance between two consecutive density matrices $\rho(t)$ and $\rho(t+dt)$ is smaller than this value. The default is tol=1e-3.

For smaller system sizes finding eigenvectors of super-operators is the preferred method:

A simple example, that can also be solved analytically is to compute the steady-state photon number of a cavity driven with an amplitude $\eta$ and damping rate $2\kappa$. The analytical solution is $n=\eta^2/\kappa^2$, while in the toolbox we may calculate this with:

b = FockBasis(10)
κ = 1.0
η = 0.5κ
H = η*(destroy(b) + create(b))
J = [sqrt(2κ)*destroy(b)]
println("n_eig = ", real(expect(number(b), ρ_eig)))
n_master = 0.24919267826355287
n_eig = 0.24999999999945127
Note, that steadystate.master returns a list of times and density matrices containing the initial time/state and the time/state when steady state is reached. Due to numerical errors the results are not exactly $0.25$. Decreasing the value of tol below its default 1e-3 improves the accuracy for the time-evolution approach.