diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..e38ba30
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,37 @@
+# Python
+__pycache__/
+*.py[cod]
+*$py.class
+*.so
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# Testing
+.coverage
+.pytest_cache/
+htmlcov/
+
+# IDEs
+.vscode/
+.idea/
+*.swp
+*.swo
+*~
+
+# OS
+.DS_Store
+Thumbs.db
diff --git a/dyson/representations/spectral.py b/dyson/representations/spectral.py
index 881561c..4092687 100644
--- a/dyson/representations/spectral.py
+++ b/dyson/representations/spectral.py
@@ -285,7 +285,7 @@ def combine(self, *args: Spectral, chempot: float | None = None) -> Spectral:
         right = util.rotate_subspace(right, orth.T.conj())
 
         # Construct the eigenvectors
-        eigvecs = util.project_eigenvectors(
+        eigvecs = util.rotate_to_physical_auxiliary_basis(
             None,
             left,
             right if not hermitian else None,
diff --git a/dyson/solvers/static/davidson.py b/dyson/solvers/static/davidson.py
index 5ae1343..e34e097 100644
--- a/dyson/solvers/static/davidson.py
+++ b/dyson/solvers/static/davidson.py
@@ -299,7 +299,7 @@ def _callback(env: dict[str, Any]) -> None:
         converged = converged[mask]
 
         # Get the full map onto physical + auxiliary and rotate the eigenvectors
-        eigvecs = util.project_eigenvectors(
+        eigvecs = util.rotate_to_physical_auxiliary_basis(
             eigvecs, self.bra, self.ket if not self.hermitian else None
         )
 
diff --git a/dyson/solvers/static/exact.py b/dyson/solvers/static/exact.py
index 3a2c091..d8cf709 100644
--- a/dyson/solvers/static/exact.py
+++ b/dyson/solvers/static/exact.py
@@ -189,7 +189,7 @@ def kernel(self) -> Spectral:
             eigvecs = np.array([left, right])
 
         # Get the full map onto physical + auxiliary and rotate the eigenvectors
-        eigvecs = util.project_eigenvectors(
+        eigvecs = util.rotate_to_physical_auxiliary_basis(
             eigvecs, self.bra, self.ket if not self.hermitian else None
         )
 
diff --git a/dyson/util/__init__.py b/dyson/util/__init__.py
index 13e96ac..9edec09 100644
--- a/dyson/util/__init__.py
+++ b/dyson/util/__init__.py
@@ -28,7 +28,7 @@
     as_trace,
     unit_vector,
     null_space_basis,
-    project_eigenvectors,
+    rotate_to_physical_auxiliary_basis,
     concatenate_paired_vectors,
     unpack_vectors,
     block_diag,
diff --git a/dyson/util/linalg.py b/dyson/util/linalg.py
index f9a214f..c9fd268 100644
--- a/dyson/util/linalg.py
+++ b/dyson/util/linalg.py
@@ -300,25 +300,44 @@ def null_space_basis(
 
 
 @cache_by_id
-def project_eigenvectors(
+def rotate_to_physical_auxiliary_basis(
     eigvecs: Array | None,
     bra: Array,
     ket: Array | None = None,
 ) -> Array:
-    """Project eigenvectors onto the physical plus auxiliary space.
+    """Rotate eigenvectors to a biorthonormal physical plus auxiliary basis.
+
+    This function constructs a complete biorthonormal basis that spans both the physical
+    space (defined by the provided bra/ket Dyson orbitals) and its orthogonal complement
+    (the auxiliary space). It then transforms the input eigenvectors into this basis.
+
+    The procedure involves:
+    1. Computing the null space of the projector formed from bra and ket vectors to obtain
+       the auxiliary space basis vectors
+    2. For Hermitian systems: Concatenating physical (bra) and auxiliary basis vectors
+    3. For non-Hermitian systems: Constructing a biorthonormal basis by orthogonalizing
+       the physical space and then biorthonormalizing the combined physical+auxiliary space
+    4. Rotating the input eigenvectors (or identity if None) into this complete basis
 
     Args:
-        eigvecs: Eigenvectors to be projected. If ``None``, assume identity.
-        bra: Bra state vector mapping the supermatrix to the physical space.
-        ket: Ket state vector mapping the supermatrix to the physical space. If ``None``, use the
-            same vectors as ``bra``.
+        eigvecs: Eigenvectors to be rotated into the physical+auxiliary basis. If ``None``,
+            returns the basis transformation itself (identity eigenvectors).
+        bra: Bra state vectors (Dyson orbitals) defining the physical space.
+            Shape: (nphys, nstates).
+        ket: Ket state vectors (Dyson orbitals) defining the physical space for non-Hermitian
+            systems. If ``None``, assumes a Hermitian system where ket = bra.
+            Shape: (nphys, nstates).
 
     Returns:
-        Projected eigenvectors.
+        Rotated eigenvectors in the physical+auxiliary basis. For Hermitian systems, shape is
+        (nphys+naux, nstates). For non-Hermitian systems, shape is (2, nphys+naux, nstates)
+        with left and right eigenvectors.
 
     Notes:
-        The physical space is defined by the ``bra`` and ``ket`` vectors, while the auxiliary part
-        is defined by the null space of the projector formed by the outer product of these vectors.
+        The physical space is defined by the ``bra`` and ``ket`` vectors (Dyson orbitals),
+        while the auxiliary space is the null space orthogonal to the projector formed by
+        the outer product of these vectors. This construction ensures the full supermatrix
+        space is spanned by a biorthonormal basis.
     """
     hermitian = ket is None
     nphys = bra.shape[0]