Fix steady-state errors

When returning the system state two errors crept into v5 that were not in v47;  the ordering of the state was not specified as Fortran ordering, and a conjugate was used instead of an adjoint in ensuring hermiticity of states.   Both are fixed, and a small test is added.

qutip/solver/heom/bofin_solvers.py

@@ -965,9 +965,10 @@ def steady_state(
         else:
             L = L.tocsc()
             solution = spsolve(L, b_mat)
-
-        data = _data.Dense(solution[:n ** 2].reshape((n, n)))
-        data = _data.mul(_data.add(data, data.conj()), 0.5)
+
+        data = _data.Dense(solution[:n ** 2].reshape((n, n), order = 'F'))                
+        data = _data.mul(_data.add(data, data.adjoint()), 0.5)
+
         steady_state = Qobj(data, dims=self._sys_dims)
 
         solution = solution.reshape((self._n_ados, n, n))

qutip/tests/solver/heom/test_bofin_solvers.py

@@ -8,8 +8,8 @@
 from scipy.integrate import quad
 
 from qutip import (
-    basis, destroy, expect, liouvillian, qeye, sigmax, sigmaz,
-    tensor, Qobj, QobjEvo
+    basis, destroy, expect, liouvillian, qeye, sigmax, sigmaz, sigmay,
+    tensor, Qobj, QobjEvo, fidelity
 )
 from qutip.core import data as _data
 from qutip.solver.heom.bofin_baths import (
@@ -740,6 +740,26 @@ def test_pure_dephasing_model_bosonic_bath(
             assert rho_final == ado_state.extract(0)
         else:
             assert_raises_steady_state_time_dependent(hsolver)
+
+
+    def test_steady_state_bosonic_bath(
+        self, atol=1e-3
+    ):
+        H_sys =  0.25 * sigmaz() + 0.5 * sigmay()
+
+        bath = DrudeLorentzBath(sigmaz(), lam=0.025, gamma=0.05, T=1/0.95, Nk=2)
+        options = {"nsteps": 15000, "store_states": True}
+        hsolver = HEOMSolver(H_sys, bath, 5, options=options)
+
+        tlist = np.linspace(0, 500, 21)
+        rho0 = basis(2, 0) * basis(2, 0).dag()
+
+        result = hsolver.run(rho0, tlist)
+        rho_final, ado_state = hsolver.steady_state()
+        fid = fidelity(rho_final, result.states[-1])
+        np.testing.assert_allclose(fid, 1.0, atol=atol)
+
+
 
     @pytest.mark.parametrize(['terminator'], [
         pytest.param(True, id="terminator"),