Repository: qutip/qutip (python)

Issue:
Dimension mismatch in solving energy restricted eigenstates
### Bug Description

Calling `Qobj.eigenstates()` of an operator created with energy restricted functionality raises an error. I believe this is due to the dimension refactoring of v5. The dimension of the solved eigenvectors are set by this list representation from `Qobj.dims` instead of the underlying `Dimension`/`Space` instances. Specifically here, the list is used instead of the `EnrSpace` instance, and thus information about the excitation limit is lost, raising a dimension/shape mismatch error in the Qobj constructor.

The simplest solution, inspired by another ENR issue #2387, would be to use `_dims` instead of `dims` in lines  1509-1512 (master at the time of writing) of `qobj.py`:
```
        if self.type == 'super':
            new_dims = [self._dims[0], [1]]
        else:
            new_dims = [self._dims[0], [1]*len(self.dims[0])]
```
I quickly checked that the superoperators work as well, but I am not sure how general the approach is. 

### Code to Reproduce the Bug

```shell
# Hamiltonian of two uncoupled qubits (harmonic oscillators limited to 2 levels) with energies e1 and e2, and only one excitation allowed in the system.

import qutip as qt
dims = [2, 2]
excitations = 1
e1, e2 = 1, 2
a1, a2 = qt.enr_destroy(dims, excitations)
H = e1*a1.dag()*a1 + e2*a2.dag()*a2
H.eigenstates()
```


### Code Output

```shell
ValueError                                Traceback (most recent call last)
Cell In[3], line 6
      4 a1, a2 = qt.enr_destroy(dims, excitations)
      5 H = e1*a1.dag()*a1 + e2*a2.dag()*a2
----> 6 H.eigenstates()

File ~\anaconda3\envs\bare_qt5\Lib\site-packages\qutip\core\qobj.py:1514, in Qobj.eigenstates(self, sparse, sort, eigvals, tol, maxiter, phase_fix)
   1512     new_dims = [self.dims[0], [1]*len(self.dims[0])]
   1513 ekets = np.empty((evecs.shape[1],), dtype=object)
-> 1514 ekets[:] = [Qobj(vec, dims=new_dims, copy=False)
   1515             for vec in _data.split_columns(evecs, False)]
   1516 norms = np.array([ket.norm() for ket in ekets])
   1517 if phase_fix is None:

File ~\anaconda3\envs\bare_qt5\Lib\site-packages\qutip\core\qobj.py:279, in Qobj.__init__(self, arg, dims, copy, superrep, isherm, isunitary)
    277 self._isherm = isherm
    278 self._isunitary = isunitary
--> 279 self._initialize_data(arg, dims, copy)
    281 if superrep is not None:
    282     self.superrep = superrep

File ~\anaconda3\envs\bare_qt5\Lib\site-packages\qutip\core\qobj.py:265, in Qobj._initialize_data(self, arg, dims, copy)
    261     self._dims = Dimensions(
    262         dims or [[self._data.shape[0]], [self._data.shape[1]]]
    263     )
    264 if self._dims.shape != self._data.shape:
--> 265     raise ValueError('Provided dimensions do not match the data: ' +
    266                      f"{self._dims.shape} vs {self._data.shape}")

ValueError: Provided dimensions do not match the data: (np.int64(4), np.int64(1)) vs (3, 1)
```


### Expected Behaviour

I expect v5 to work like v4.7.5, where the same code gives the following output:
``` 
(array([0., 1., 2.]),
 array([Quantum object: dims = [[2, 2], [1, 1]], shape = (3, 1), type = ket
        Qobj data =
        [[1.]
         [0.]
         [0.]]                                                             ,
        Quantum object: dims = [[2, 2], [1, 1]], shape = (3, 1), type = ket
        Qobj data =
        [[0.]
         [0.]
         [1.]]                                                             ,
        Quantum object: dims = [[2, 2], [1, 1]], shape = (3, 1), type = ket
        Qobj data =
        [[0.]
         [1.]
         [0.]]                                                             ],
       dtype=object))
```


### Your Environment

```shell
QuTiP Version:      5.1.0
Numpy Version:      2.2.1
Scipy Version:      1.15.0
Cython Version:     None
Matplotlib Version: 3.10.0
Python Version:     3.12.8
Number of CPUs:     8
BLAS Info:          Generic
INTEL MKL Ext:      None
Platform Info:      Windows (AMD64)
```


### Additional Context

Contributing to the Qutip community is new territory for me, so please bare with me as learn the ropes. I'm happy to create a PR with additional testing for this issue, if the suggested solution seems reasonable.

Pull request:
ENR dimensions work with eigenstates
**Description**
Use the class dims instead of the list dims in `eigenstates` so it support ENR.

**Related issues or PRs**
fix #2595


Patch:
diff --git a/doc/changes/2595.bugfix b/doc/changes/2595.bugfix
new file mode 100644
index 000000000..73ce1215a
--- /dev/null
+++ b/doc/changes/2595.bugfix
@@ -0,0 +1,1 @@
+ENR dimensions work with eigenstates
diff --git a/qutip/core/qobj.py b/qutip/core/qobj.py
index 811477120..0692cc003 100644
--- a/qutip/core/qobj.py
+++ b/qutip/core/qobj.py
@@ -1507,9 +1507,9 @@ class Qobj:
                                       sort=sort, eigvals=eigvals)
 
         if self.type == 'super':
-            new_dims = [self.dims[0], [1]]
+            new_dims = [self._dims[0], [1]]
         else:
-            new_dims = [self.dims[0], [1]*len(self.dims[0])]
+            new_dims = [self._dims[0], [1]*len(self.dims[0])]
         ekets = np.empty((evecs.shape[1],), dtype=object)
         ekets[:] = [Qobj(vec, dims=new_dims, copy=False)
                     for vec in _data.split_columns(evecs, False)]


Test patch:
diff --git a/qutip/tests/test_enr_state_operator.py b/qutip/tests/test_enr_state_operator.py
index 9257e3266..3919151fc 100644
--- a/qutip/tests/test_enr_state_operator.py
+++ b/qutip/tests/test_enr_state_operator.py
@@ -189,3 +189,11 @@ def test_steadystate_ENR():
 
     np.testing.assert_allclose(exp_sz_JC,
                                exp_sz_enr, atol=1e-2)
+
+
+def test_eigenstates_ENR():
+    a1, a2 = qutip.enr_destroy([2, 2], 1)
+    H = a1.dag() * a1 + a2.dag() * a2
+    vals, vecs = H.eigenstates()
+    for val, vec in zip(vals, vecs):
+        assert val * vec == H @ vec