  Symmetries in an open quantum system lead to degenerated Liouvillian that
physically implies the existence of multiple steady states. In such cases,
obtaining the initial condition independent stead states is highly nontrivial
since any linear combination of the \emph{true} asymptotic states, which may
not necessarily be a density matrix, is also a valid asymptote for the
Liouvillian. Thus, in this work we consider different approaches to obtain the
\emph{true} steady states of a degenerated Liouvillian. In the ideal scenario,
when the open system symmetry operators are known we show how these can be used
to obtain the invariant subspaces of the Liouvillian and hence the steady
states. We then discuss two other approaches that do not require any knowledge
of the symmetry operators. These could be a powerful tool to deal with quantum
many-body complex open systems. The first approach which is based on
Gramm-Schmidt orthonormalization of density matrices allows us to obtain
\emph{all} the steady states, whereas the second one based on large deviations
allows us to obtain the non-degenerated maximum and minimum current-carrying
states. We discuss our method with the help of an open para-Benzene ring and
examine interesting scenarios such as the dynamical restoration of Hamiltonian
symmetries in the long-time limit and apply the method to study the
eigenspacing statistics of the nonequilibrium steady state.
