diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
index 36070231..d0f64e72 100644
--- a/.github/workflows/ci.yaml
+++ b/.github/workflows/ci.yaml
@@ -161,11 +161,8 @@ jobs:
           python-version: ${{ matrix.version }}
           auto-activate-base: false
       - name: Install Psi4 from c-f
-        # TODO: remove libxc-c from the install list once
-        # https://github.com/psi4/psi4/issues/3474
-        # is fixed
         run: |
-          mamba install psi4 pyddx "libxc-c=7.0.0" -c conda-forge -c conda-forge/label/libint_dev
+          mamba install psi4 pyddx -c conda-forge -c conda-forge/label/libint_dev
           psi4 --version
           python -c "import psi4"
       - name: Install system dependencies on macOS
diff --git a/adcc/backends/EriBuilder.py b/adcc/backends/EriBuilder.py
index 7e5e7eb2..258e8d96 100644
--- a/adcc/backends/EriBuilder.py
+++ b/adcc/backends/EriBuilder.py
@@ -20,22 +20,47 @@
 ## along with adcc. If not, see <http://www.gnu.org/licenses/>.
 ##
 ## ---------------------------------------------------------------------
+from dataclasses import dataclass
 from itertools import product
-from collections import namedtuple
+from typing import Literal, TypeAlias, TypeGuard
+import numpy as np
 
+IntSlice: TypeAlias = "slice[int, int, int]"
+IntSlice4D = tuple[IntSlice, IntSlice, IntSlice, IntSlice]
+Block = Literal["O", "V"]
+Block4D = tuple[Block, Block, Block, Block]
+Spin = Literal["a", "b"]
+Spin4D = tuple[Spin, Spin, Spin, Spin]
+Array4D = np.ndarray[tuple[int, int, int, int], np.dtype[np.float64]]
 
-def range_in(inner, full):
-    if inner.start is None:
-        inner = slice(0, inner.stop, 1)
-    if full.start is None:
-        full = slice(0, full.stop, 1)
-    return all(r in range(full.start, full.stop)
-               for r in range(inner.start, inner.stop))
+
+def is_int_slice(slice: slice) -> TypeGuard[IntSlice]:
+    return (
+        isinstance(slice.start, int)
+        and isinstance(slice.stop, int)
+        and isinstance(slice.step, int)
+    )
 
 
-# Helper namedtuple for slices of spin blocks
-SpinBlockSlice = namedtuple('SpinBlockSlice',
-                            ['block', 'spin', 'fromslice', 'toslice'])
+@dataclass(frozen=True, slots=True)
+class SpinBlockSlice:
+    block: Block
+    spin: Spin
+    fromslice: IntSlice
+    toslice: IntSlice
+
+
+@dataclass(frozen=True, slots=True)
+class SpinBlockSlice4D:
+    block: Block4D
+    spin: Spin4D
+    fromslice: IntSlice4D
+    toslice: IntSlice4D
+
+
+def range_in(inner: IntSlice, full: IntSlice) -> bool:
+    return all(r in range(full.start, full.stop)
+               for r in range(inner.start, inner.stop))
 
 
 class EriBuilder:
@@ -48,49 +73,56 @@ class EriBuilder:
           Gets passed the block as a string like 'OOVV' and the spin block as
           as string like 'abab'.
     """
-    def __init__(self, n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted):
-        self.n_orbs = n_orbs
-        self.n_orbs_alpha = n_orbs_alpha
-        self.n_alpha = n_alpha
-        self.n_beta = n_beta
-        self.eri_cache = {}
-        self.restricted = restricted
-        self.block2slice = {
-            "oa": slice(0, self.n_alpha, 1),
-            "va": slice(self.n_alpha, self.n_orbs_alpha, 1),
-            "ob": slice(self.n_orbs_alpha, self.n_orbs_alpha + self.n_beta, 1),
-            "vb": slice(self.n_orbs_alpha + self.n_beta, self.n_orbs, 1),
+    def __init__(self, n_orbs: int, n_orbs_alpha: int, n_alpha: int, n_beta: int,
+                 restricted: bool):
+        self.n_orbs: int = n_orbs
+        self.n_orbs_alpha: int = n_orbs_alpha
+        self.n_alpha: int = n_alpha
+        self.n_beta: int = n_beta
+        self.eri_cache: dict[str, Array4D] = {}
+        self.restricted: bool = restricted
+        self.block2slice: dict[tuple[Block, Spin], IntSlice] = {
+            ("O", "a"): slice(0, self.n_alpha, 1),
+            ("V", "a"): slice(self.n_alpha, self.n_orbs_alpha, 1),
+            ("O", "b"): slice(self.n_orbs_alpha,
+                              self.n_orbs_alpha + self.n_beta, 1),
+            ("V", "b"): slice(self.n_orbs_alpha + self.n_beta, self.n_orbs, 1),
         }
 
-    def compute_mo_eri(self, blocks, spins):
+    def compute_mo_eri(self, blocks: Block4D, spins: Spin4D) -> Array4D:
         """
         Compute block of the ERI tensor in chemists' indexing
         """
         raise NotImplementedError("Implement compute_mo_eri")
 
-    def split_4d_slice(self, slices):
+    def split_4d_slice(
+        self, slices: tuple[slice, slice, slice, slice]
+    ) -> list[SpinBlockSlice4D]:
         """
         Split tuple of four slices into the block spin slices
         and their mapping to where elements are to be placed
         """
-        return [SpinBlockSlice(tpl[0][0] + tpl[1][0] + tpl[2][0] + tpl[3][0],
-                               tpl[0][1] + tpl[1][1] + tpl[2][1] + tpl[3][1],
-                               (tpl[0][2], tpl[1][2], tpl[2][2], tpl[3][2]),
-                               (tpl[0][3], tpl[1][3], tpl[2][3], tpl[3][3]))
-                for tpl in product(*(self.split_1d_slice(sl) for sl in slices))]
+        splitted = (self.split_1d_slice(sl) for sl in slices)
+        return [SpinBlockSlice4D(
+            (sl1.block, sl2.block, sl3.block, sl4.block),
+            (sl1.spin, sl2.spin, sl3.spin, sl4.spin),
+            (sl1.fromslice, sl2.fromslice, sl3.fromslice, sl4.fromslice),
+            (sl1.toslice, sl2.toslice, sl3.toslice, sl4.toslice)
+        ) for sl1, sl2, sl3, sl4 in product(*splitted)]
 
-    def split_1d_slice(self, sl):
+    def split_1d_slice(self, sl: slice) -> list[SpinBlockSlice]:
         """
         Split slice into block-slices or multiple block-slices
         """
         if sl.start is None:
-            sl = slice(0, sl.stop, 1)
+            sl = slice(0, sl.stop, sl.step)
         if sl.step is None:
             sl = slice(sl.start, sl.stop, 1)
-
-        ret = []
+        assert is_int_slice(sl)
+        ret: list[SpinBlockSlice] = []
         for (block, bslice) in self.block2slice.items():
-            fromslice = toslice = None
+            fromslice: tuple[int, int] | None = None
+            toslice: tuple[int, int] | None = None
             if range_in(sl, bslice):
                 fromslice = (sl.start - bslice.start, sl.stop - bslice.start)
                 toslice = (0, sl.stop - sl.start)
@@ -107,20 +139,32 @@ def split_1d_slice(self, sl):
                 toslice = (bslice.start - sl.start, sl.stop - sl.start)
             if fromslice is None or toslice is None:
                 continue   # Not found
-            ret.append(SpinBlockSlice(block[0].upper(), block[1],
-                                      slice(*fromslice), slice(*toslice)))
+            ret.append(SpinBlockSlice(
+                block[0], block[1], slice(*fromslice, 1), slice(*toslice, 1)
+            ))
         assert len(ret) > 0
         return ret
 
-    def fill_slice_symm(self, slices, out):
+    def fill_slice_symm(
+        self, slices: tuple[slice, slice, slice, slice], out: Array4D
+    ) -> None:
+        non_zero_spin_blocks: list[Spin4D] = [  # chemist notation
+            ("a", "a", "a", "a"),
+            ("a", "a", "b", "b"),
+            ("b", "b", "a", "a"),
+            ("b", "b", "b", "b"),
+        ]
         for sbslices in self.split_4d_slice(slices):
-            blocks, spins, fromslices, toslices = sbslices
-            if spins not in ["aaaa", "aabb", "bbaa", "bbbb"]:
+            blocks: Block4D = sbslices.block
+            spins: Spin4D = sbslices.spin
+            fromslices: IntSlice4D = sbslices.fromslice
+            toslices: IntSlice4D = sbslices.toslice
+            if spins not in non_zero_spin_blocks:
                 out[toslices] = 0  # Zero by symmetry
                 continue
             if self.restricted:
                 # For restricted spins in chem eri do not matter
-                spins = "aaaa"
+                spins = ("a", "a", "a", "a")
 
             cache_key = "".join(blocks) + "".join(spins)
             if cache_key in self.eri_cache:
@@ -131,5 +175,5 @@ def fill_slice_symm(self, slices, out):
 
             out[toslices] = eri[fromslices]
 
-    def flush_cache(self):
+    def flush_cache(self) -> None:
         self.eri_cache = {}
diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index cb92a348..76974537 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -20,16 +20,32 @@
 ## along with adcc. If not, see <http://www.gnu.org/licenses/>.
 ##
 ## ---------------------------------------------------------------------
+from typing import Literal
 import numpy as np
 
-from libadcc import HartreeFockProvider
-
 import psi4
 
-from .EriBuilder import EriBuilder
+import libadcc
+
+from .EriBuilder import EriBuilder, Block4D, Spin4D
 from ..exceptions import InvalidReference
 from ..ElectronicStates import EnergyCorrection
-from ..OneParticleOperator import OneParticleOperator
+
+# Some type defs for the interface
+Array1D = np.ndarray[tuple[int], np.dtype[np.float64]]
+Array2D = np.ndarray[tuple[int, int], np.dtype[np.float64]]
+Array4D = np.ndarray[tuple[int, int, int, int], np.dtype[np.float64]]
+DipoleLike = tuple[Array2D, Array2D, Array2D]
+# Once we drop python 3.10 we can write
+# QuadrupoleLike = tuple[*DipoleLike, *DipoleLike, *DipoleLike]
+QuadrupoleLike = tuple[
+    Array2D, Array2D, Array2D,
+    Array2D, Array2D, Array2D,
+    Array2D, Array2D, Array2D,
+]
+Coordinate = tuple[float, float, float]
+Environment = Literal["pe", "pcm"]
+EnvironmentImplementation = Literal["cppe", "ddx", "pcmsolver"]
 
 
 class Psi4OperatorIntegralProvider:
@@ -39,21 +55,23 @@ class Psi4OperatorIntegralProvider:
         "pe_induction_elec", "pcm_potential_elec"
     )
 
-    def __init__(self, wfn):
-        self.wfn = wfn
-        self.backend = "psi4"
-        self.mints = psi4.core.MintsHelper(self.wfn)
+    def __init__(self, wfn: psi4.core.HF):
+        self.wfn: psi4.core.HF = wfn
+        self.backend: str = "psi4"
+        self.mints: psi4.core.MintsHelper = psi4.core.MintsHelper(self.wfn)
 
     @property
-    def overlap(self) -> np.ndarray:
+    def overlap(self) -> Array2D:
         return np.asarray(self.mints.ao_overlap())
 
     @property
-    def electric_dipole(self) -> tuple[np.ndarray, ...]:
+    def electric_dipole(self) -> DipoleLike:
         """-sum_i r_i"""
-        return tuple(np.asarray(comp) for comp in self.mints.ao_dipole())
+        x, y, z = self.mints.ao_dipole()  # list
+        return np.asarray(x), np.asarray(y), np.asarray(z)
 
-    def magnetic_dipole(self, gauge_origin="origin") -> tuple[np.ndarray, ...]:
+    def magnetic_dipole(self,
+                        gauge_origin: Coordinate | str = "origin") -> DipoleLike:
         """
         The imaginary part of the integral is returned.
         -0.5 * sum_i r_i x p_i
@@ -65,20 +83,21 @@ def magnetic_dipole(self, gauge_origin="origin") -> tuple[np.ndarray, ...]:
                 " gauge origin for the magnetic dipole operator. "
                 f"{gauge_origin} is not valid."
             )
-        return tuple(
-            0.5 * np.asarray(comp)
-            for comp in self.mints.ao_angular_momentum()
-        )
+        x, y, z = self.mints.ao_angular_momentum()  # list
+        return 0.5 * np.asarray(x), 0.5 * np.asarray(y), 0.5 * np.asarray(z)
 
     @property
-    def electric_dipole_velocity(self) -> tuple[np.ndarray, ...]:
+    def electric_dipole_velocity(self) -> DipoleLike:
         """
         The imaginary part of the integral is returned.
         -sum_i p_i
         """
-        return tuple(-1.0 * np.asarray(comp) for comp in self.mints.ao_nabla())
+        x, y, z = self.mints.ao_nabla()  # list
+        return -1.0 * np.asarray(x), -1.0 * np.asarray(y), -1.0 * np.asarray(z)
 
-    def electric_quadrupole(self, gauge_origin="origin") -> tuple[np.ndarray, ...]:
+    def electric_quadrupole(
+        self, gauge_origin: Coordinate | str = "origin"
+    ) -> QuadrupoleLike:
         """-sum_i r_{i, alpha} r_{i, beta}"""
         # TODO: Gauge origin?
         if gauge_origin != (0.0, 0.0, 0.0) and gauge_origin != "origin":
@@ -92,8 +111,9 @@ def electric_quadrupole(self, gauge_origin="origin") -> tuple[np.ndarray, ...]:
         assert len(u) == 6
         return (u[0], u[1], u[2], u[1], u[3], u[4], u[2], u[4], u[5])
 
-    def electric_quadrupole_traceless(self, gauge_origin="origin"
-                                      ) -> tuple[np.ndarray, ...]:
+    def electric_quadrupole_traceless(
+        self, gauge_origin: Coordinate | str = "origin"
+    ) -> QuadrupoleLike:
         """
         -0.5 * sum_i (3 * r_{i, alpha} r_{i, beta}
         - delta_{alpha, beta} r_{i}^2)
@@ -110,39 +130,41 @@ def electric_quadrupole_traceless(self, gauge_origin="origin"
         assert len(u) == 6
         return (u[0], u[1], u[2], u[1], u[3], u[4], u[2], u[4], u[5])
 
-    def pe_induction_elec(self, dm: OneParticleOperator) -> psi4.core.Matrix:
+    def pe_induction_elec(self, dm: libadcc.Tensor) -> Array2D:
         if not hasattr(self.wfn, "pe_state"):
             raise RuntimeError("Can not compute the PE electronic induction "
                                "operator using the given psi4 object.")
 
-        return self.wfn.pe_state.get_pe_contribution(
+        return np.asarray(self.wfn.pe_state.get_pe_contribution(
             psi4.core.Matrix.from_array(dm.to_ndarray()),
             elec_only=True
-        )[1]
+        )[1])
 
-    def pcm_potential_elec(self, dm: OneParticleOperator) -> psi4.core.Matrix:
+    def pcm_potential_elec(self, dm: libadcc.Tensor) -> Array2D:
         if hasattr(self.wfn, "ddx"):
-            return self.wfn.ddx.get_solvation_contributions(
+            return np.asarray(self.wfn.ddx.get_solvation_contributions(
                 psi4.core.Matrix.from_array(dm.to_ndarray()),
                 elec_only=True, nonequilibrium=True
-            )[1]
+            )[1])
         elif self.wfn.PCM_enabled():
-            return psi4.core.PCM.compute_V(
-                self.wfn.get_PCM(),
+            return np.asarray(self.wfn.get_PCM().compute_V(
                 psi4.core.Matrix.from_array(dm.to_ndarray())
-            )
+            ))
         raise RuntimeError("Can not compute the electronic PCM potential "
                            "operator using the given psi4 object.")
 
 
 class Psi4EriBuilder(EriBuilder):
-    def __init__(self, wfn, n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted):
-        self.wfn = wfn
-        self.mints = psi4.core.MintsHelper(self.wfn)
+    def __init__(self, wfn: psi4.core.HF, n_orbs: int, n_orbs_alpha: int,
+                 n_alpha: int, n_beta: int, restricted: bool):
+        self.wfn: psi4.core.HF = wfn
+        self.mints: psi4.core.MintsHelper = psi4.core.MintsHelper(self.wfn)
         super().__init__(n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted)
 
     @property
-    def coefficients(self):
+    def coefficients(
+        self
+    ) -> dict[Literal["Oa", "Ob", "Va", "Vb"], psi4.core.Matrix]:
         return {
             "Oa": self.wfn.Ca_subset("AO", "OCC"),
             "Ob": self.wfn.Cb_subset("AO", "OCC"),
@@ -150,28 +172,31 @@ def coefficients(self):
             "Vb": self.wfn.Cb_subset("AO", "VIR"),
         }
 
-    def compute_mo_eri(self, blocks, spins):
+    def compute_mo_eri(self, blocks: Block4D, spins: Spin4D) -> Array4D:
         coeffs = tuple(self.coefficients[blocks[i] + spins[i]] for i in range(4))
         return np.asarray(self.mints.mo_eri(*coeffs))
 
 
-class Psi4HFProvider(HartreeFockProvider):
+class Psi4HFProvider(libadcc.HartreeFockProvider):
     """
         This implementation is only valid
         if no orbital reordering is required.
     """
-    def __init__(self, wfn):
+    def __init__(self, wfn: psi4.core.HF):
         # Do not forget the next line,
         # otherwise weird errors result
         super().__init__()
-        self.wfn = wfn
-        self.eri_builder = Psi4EriBuilder(self.wfn, self.n_orbs, self.wfn.nmo(),
-                                          wfn.nalpha(), wfn.nbeta(),
-                                          self.restricted)
-        self.operator_integral_provider = Psi4OperatorIntegralProvider(self.wfn)
-
-        self.environment = None
-        self.environment_implementation = None
+        self.wfn: psi4.core.HF = wfn
+        self.eri_builder: Psi4EriBuilder = Psi4EriBuilder(
+            self.wfn, self.n_orbs, self.wfn.nmo(),
+            wfn.nalpha(), wfn.nbeta(), self.restricted
+        )
+        self.operator_integral_provider: Psi4OperatorIntegralProvider = (
+            Psi4OperatorIntegralProvider(self.wfn)
+        )
+
+        self.environment: Environment | None = None
+        self.environment_implementation: EnvironmentImplementation | None = None
         if hasattr(self.wfn, "pe_state"):
             self.environment = "pe"
             self.environment_implementation = "cppe"
@@ -182,13 +207,13 @@ def __init__(self, wfn):
             self.environment = "pcm"
             self.environment_implementation = "pcmsolver"
 
-    def pe_energy(self, dm, elec_only=True):
+    def pe_energy(self, dm: libadcc.Tensor, elec_only: bool = True) -> float:
         density_psi = psi4.core.Matrix.from_array(dm.to_ndarray())
         e_pe, _ = self.wfn.pe_state.get_pe_contribution(density_psi,
                                                         elec_only=elec_only)
         return e_pe
 
-    def pcm_energy(self, dm, elec_only=True):
+    def pcm_energy(self, dm: libadcc.Tensor, elec_only: bool = True) -> float:
         psi_dm = psi4.core.Matrix.from_array(dm.to_ndarray())
         # computes the Fock matrix contribution.
         # By contraction with the tdm, the electronic energy contribution is
@@ -199,10 +224,13 @@ def pcm_energy(self, dm, elec_only=True):
             )[1]
         elif self.environment_implementation == "pcmsolver":
             V_pcm = psi4.core.PCM.compute_V(self.wfn.get_PCM(), psi_dm)
-        return np.einsum("uv,uv->", dm.to_ndarray(), V_pcm.to_array())
+        else:
+            raise ValueError("Invalid environment implementation: "
+                             f"{self.environment_implementation}")
+        return float(np.einsum("uv,uv->", dm.to_ndarray(), V_pcm.to_array()))
 
     @property
-    def excitation_energy_corrections(self):
+    def excitation_energy_corrections(self) -> dict[str, EnergyCorrection]:
         ret = []
         if self.environment == "pe":
             ptlr = EnergyCorrection(
@@ -225,35 +253,38 @@ def excitation_energy_corrections(self):
             ret.extend([ptlr])
         return {ec.name: ec for ec in ret}
 
-    def get_backend(self):
+    def get_backend(self) -> str:
         return "psi4"
 
-    def get_conv_tol(self):
+    def get_conv_tol(self) -> float:
         conv_tol = psi4.core.get_option("SCF", "E_CONVERGENCE")
         # RMS value of the orbital gradient
         conv_tol_grad = psi4.core.get_option("SCF", "D_CONVERGENCE")
         threshold = max(conv_tol, conv_tol_grad)
         return threshold
 
-    def get_restricted(self):
+    def get_restricted(self) -> bool:
         return isinstance(self.wfn, (psi4.core.RHF, psi4.core.ROHF))
 
-    def get_energy_scf(self):
+    def get_energy_scf(self) -> float:
         return self.wfn.energy()
 
-    def get_nuclear_repulsion_energy(self):
+    def get_nuclear_repulsion_energy(self) -> float:
         return self.wfn.molecule().nuclear_repulsion_energy()
 
-    def get_spin_multiplicity(self):
+    def get_spin_multiplicity(self) -> int:
         return self.wfn.molecule().multiplicity()
 
-    def get_n_orbs_alpha(self):
+    def get_n_orbs_alpha(self) -> int:
         return self.wfn.nmo()
 
-    def get_n_bas(self):
+    def get_n_bas(self) -> int:
         return self.wfn.basisset().nbf()
 
-    def get_nuclear_multipole(self, order, gauge_origin=(0, 0, 0)):
+    def get_nuclear_multipole(
+        self, order: int,
+        gauge_origin: Coordinate = (0.0, 0.0, 0.0)
+    ) -> Array1D:
         molecule = self.wfn.molecule()
         if order == 0:
             # The function interface needs to be a np.array on return
@@ -265,10 +296,11 @@ def get_nuclear_multipole(self, order, gauge_origin=(0, 0, 0)):
         else:
             raise NotImplementedError("get_nuclear_multipole with order > 1")
 
-    def transform_gauge_origin_to_xyz(self, gauge_origin):
+    def transform_gauge_origin_to_xyz(self, gauge_origin: str
+                                      ) -> Coordinate:
         raise NotImplementedError("transform_gauge_origin_to_xyz not implemented.")
 
-    def fill_orbcoeff_fb(self, out):
+    def fill_orbcoeff_fb(self, out: Array2D) -> None:
         mo_coeff_a = np.asarray(self.wfn.Ca())
         mo_coeff_b = np.asarray(self.wfn.Cb())
         mo_coeff = (mo_coeff_a, mo_coeff_b)
@@ -276,36 +308,38 @@ def fill_orbcoeff_fb(self, out):
             np.hstack((mo_coeff[0], mo_coeff[1]))
         )
 
-    def fill_occupation_f(self, out):
+    def fill_occupation_f(self, out: Array1D) -> None:
         out[:] = np.hstack((
             np.asarray(self.wfn.occupation_a()),
             np.asarray(self.wfn.occupation_b())
         ))
 
-    def fill_orben_f(self, out):
+    def fill_orben_f(self, out: Array1D) -> None:
         orben_a = np.asarray(self.wfn.epsilon_a())
         orben_b = np.asarray(self.wfn.epsilon_b())
         out[:] = np.hstack((orben_a, orben_b))
 
-    def fill_fock_ff(self, slices, out):
+    def fill_fock_ff(self, slices: tuple[slice, slice], out: Array2D) -> None:
         diagonal = np.empty(self.n_orbs)
         self.fill_orben_f(diagonal)
         out[:] = np.diag(diagonal)[slices]
 
-    def fill_eri_ffff(self, slices, out):
+    def fill_eri_ffff(self, slices: tuple[slice, slice, slice, slice],
+                      out: Array4D) -> None:
         self.eri_builder.fill_slice_symm(slices, out)
 
-    def fill_eri_phys_asym_ffff(self, slices, out):
+    def fill_eri_phys_asym_ffff(self, slices: tuple[slice, slice, slice, slice],
+                                out: Array4D) -> None:
         raise NotImplementedError("fill_eri_phys_asym_ffff not implemented.")
 
-    def has_eri_phys_asym_ffff(self):
+    def has_eri_phys_asym_ffff(self) -> bool:
         return False
 
-    def flush_cache(self):
+    def flush_cache(self) -> None:
         self.eri_builder.flush_cache()
 
 
-def import_scf(wfn):
+def import_scf(wfn: psi4.core.HF) -> Psi4HFProvider:
     if not isinstance(wfn, psi4.core.HF):
         raise InvalidReference(
             "Only psi4.core.HF and its subtypes are supported references in "
@@ -338,8 +372,9 @@ def import_scf(wfn):
     return provider
 
 
-def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
-           conv_tol_grad=1e-9, max_iter=150, pe_options=None):
+def run_hf(xyz: str, basis: str, charge: int = 0, multiplicity: int = 1,
+           conv_tol: float = 1e-11, conv_tol_grad: float = 1e-9,
+           max_iter: int = 150, pe_options: dict | None = None) -> psi4.core.HF:
     basissets = {
         "sto3g": "sto-3g",
         "def2tzvp": "def2-tzvp",
@@ -377,6 +412,6 @@ def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
             'soscf': 'true'
         })
 
-    _, wfn = psi4.energy('SCF', return_wfn=True, molecule=mol)
+    _, wfn = psi4.energy('SCF', return_wfn=True, molecule=mol)  # type: ignore
     psi4.core.clean()
     return wfn
diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py
index 13075826..22c52513 100644
--- a/adcc/backends/pyscf.py
+++ b/adcc/backends/pyscf.py
@@ -20,17 +20,33 @@
 ## along with adcc. If not, see <http://www.gnu.org/licenses/>.
 ##
 ## ---------------------------------------------------------------------
+from typing import cast, Literal
 import numpy as np
 
-from libadcc import HartreeFockProvider
+from pyscf import ao2mo, gto, scf
+from pyscf.solvent import ddcosmo
+
+import libadcc
 
-from .EriBuilder import EriBuilder
+from .EriBuilder import EriBuilder, Block4D, Spin4D
 from ..exceptions import InvalidReference
 from ..ElectronicStates import EnergyCorrection
-from ..OneParticleOperator import OneParticleOperator
 
-from pyscf import ao2mo, gto, scf
-from pyscf.solvent import ddcosmo
+# Some type defs for the interface
+Array1D = np.ndarray[tuple[int], np.dtype[np.float64]]
+Array2D = np.ndarray[tuple[int, int], np.dtype[np.float64]]
+Array4D = np.ndarray[tuple[int, int, int, int], np.dtype[np.float64]]
+DipoleLike = tuple[Array2D, Array2D, Array2D]
+# Once we drop python 3.10 we can write
+# QuadrupoleLike = tuple[*DipoleLike, *DipoleLike, *DipoleLike]
+QuadrupoleLike = tuple[
+    Array2D, Array2D, Array2D,
+    Array2D, Array2D, Array2D,
+    Array2D, Array2D, Array2D,
+]
+Coordinate = tuple[float, float, float]
+Environment = Literal["pe", "pcm"]
+EnvironmentImplementation = Literal["cppe", "ddcosmo"]
 
 
 class PyScfOperatorIntegralProvider:
@@ -41,22 +57,24 @@ class PyScfOperatorIntegralProvider:
         "pe_induction_elec", "pcm_potential_elec"
     )
 
-    def __init__(self, scfres):
-        self.scfres = scfres
-        self.backend = "pyscf"
+    def __init__(self, scfres: scf.hf.SCF):
+        self.scfres: scf.hf.SCF = scfres
+        self.backend: str = "pyscf"
 
     @property
-    def overlap(self) -> np.ndarray:
+    def overlap(self) -> Array2D:
         return self.scfres.mol.intor_symmetric('int1e_ovlp')
 
     @property
-    def electric_dipole(self) -> tuple[np.ndarray, ...]:
+    def electric_dipole(self) -> DipoleLike:
         """-sum_i r_i"""
         return tuple(
             -1.0 * self.scfres.mol.intor_symmetric('int1e_r', comp=3)
         )
 
-    def magnetic_dipole(self, gauge_origin="origin") -> tuple[np.ndarray, ...]:
+    def magnetic_dipole(
+        self, gauge_origin: Coordinate | str = "origin"
+    ) -> DipoleLike:
         """
         The imaginary part of the integral is returned.
         -0.5 * sum_i r_i x p_i
@@ -68,7 +86,7 @@ def magnetic_dipole(self, gauge_origin="origin") -> tuple[np.ndarray, ...]:
             )
 
     @property
-    def electric_dipole_velocity(self) -> tuple[np.ndarray, ...]:
+    def electric_dipole_velocity(self) -> DipoleLike:
         """
         The imaginary part of the integral is returned.
         -sum_i p_i
@@ -78,7 +96,9 @@ def electric_dipole_velocity(self) -> tuple[np.ndarray, ...]:
                 self.scfres.mol.intor('int1e_ipovlp', comp=3, hermi=2)
             )
 
-    def electric_quadrupole(self, gauge_origin="origin") -> tuple[np.ndarray, ...]:
+    def electric_quadrupole(
+        self, gauge_origin: Coordinate | str = "origin"
+    ) -> QuadrupoleLike:
         """-sum_i r_{i, alpha} r_{i, beta}"""
         gauge_origin = _transform_gauge_origin_to_xyz(self.scfres, gauge_origin)
         with self.scfres.mol.with_common_orig(gauge_origin):
@@ -86,8 +106,9 @@ def electric_quadrupole(self, gauge_origin="origin") -> tuple[np.ndarray, ...]:
                 -1.0 * self.scfres.mol.intor_symmetric('int1e_rr', comp=9)
             )
 
-    def electric_quadrupole_traceless(self, gauge_origin="origin"
-                                      ) -> tuple[np.ndarray, ...]:
+    def electric_quadrupole_traceless(
+        self, gauge_origin: Coordinate | str = "origin"
+    ) -> QuadrupoleLike:
         """
         -0.5 * sum_i (3 * r_{i, alpha} r_{i, beta}
         - delta_{alpha, beta} r_{i}^2)
@@ -105,8 +126,9 @@ def electric_quadrupole_traceless(self, gauge_origin="origin"
                 -1.0 * np.reshape(term, (9, r_quadr.shape[0], r_quadr.shape[0]))
             )
 
-    def electric_quadrupole_velocity(self, gauge_origin="origin"
-                                     ) -> tuple[np.ndarray, ...]:
+    def electric_quadrupole_velocity(
+        self, gauge_origin: Coordinate | str = "origin"
+    ) -> QuadrupoleLike:
         """
         The imaginary part of the integral is returned.
         -sum_i (r_{i, beta} p_{i, alpha} - i delta_{alpha, beta}
@@ -125,8 +147,9 @@ def electric_quadrupole_velocity(self, gauge_origin="origin"
                 -1.0 * np.reshape(term, (9, ovlp.shape[0], ovlp.shape[0]))
             )
 
-    def diamagnetic_magnetizability(self, gauge_origin="origin"
-                                    ) -> tuple[np.ndarray, ...]:
+    def diamagnetic_magnetizability(
+        self, gauge_origin: Coordinate | str = "origin"
+    ) -> QuadrupoleLike:
         """
         0.25 * sum_i (r_{i, alpha} r_{i, beta}
         - delta_{alpha, beta} r_{i}^2)
@@ -140,11 +163,11 @@ def diamagnetic_magnetizability(self, gauge_origin="origin"
             for i in range(3):
                 r_quadr_matrix[i][i] = r_quadr
             term = 0.25 * (r_quadr_matrix - r_r)
-            return tuple(
+            return cast(QuadrupoleLike, tuple(
                 np.reshape(term, (9, r_quadr.shape[0], r_quadr.shape[0]))
-            )
+            ))
 
-    def pe_induction_elec(self, dm: OneParticleOperator) -> np.ndarray:
+    def pe_induction_elec(self, dm: libadcc.Tensor) -> Array2D:
         try:
             self.scfres.with_solvent.cppe_state
         except AttributeError:
@@ -155,7 +178,7 @@ def pe_induction_elec(self, dm: OneParticleOperator) -> np.ndarray:
             dm.to_ndarray(), elec_only=True
         )[1]
 
-    def pcm_potential_elec(self, dm: OneParticleOperator) -> np.ndarray:
+    def pcm_potential_elec(self, dm: libadcc.Tensor) -> Array2D:
         if not hasattr(self.scfres, "with_solvent") or \
                 not isinstance(self.scfres.with_solvent, ddcosmo.DDCOSMO):
             raise RuntimeError("Can not compute the electronic PCM potential "
@@ -172,16 +195,23 @@ def pcm_potential_elec(self, dm: OneParticleOperator) -> np.ndarray:
 # TODO: refactor ERI builder to be more general
 # IntegralBuilder would be good
 class PyScfEriBuilder(EriBuilder):
-    def __init__(self, scfres, n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted):
-        self.scfres = scfres
+    def __init__(self, scfres: scf.hf.SCF, n_orbs: int, n_orbs_alpha: int,
+                 n_alpha: int, n_beta: int, restricted: bool):
+        self.scfres: scf.hf.SCF = scfres
+        self.mo_coeff: tuple[Array2D, Array2D]
         if restricted:
-            self.mo_coeff = (self.scfres.mo_coeff, self.scfres.mo_coeff)
+            self.mo_coeff = cast(
+                tuple[Array2D, Array2D],
+                (self.scfres.mo_coeff, self.scfres.mo_coeff)
+            )
         else:
-            self.mo_coeff = self.scfres.mo_coeff
+            self.mo_coeff = cast(
+                tuple[Array2D, Array2D], self.scfres.mo_coeff
+            )
         super().__init__(n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted)
 
     @property
-    def coefficients(self):
+    def coefficients(self) -> dict[Literal["Oa", "Ob", "Va", "Vb"], Array2D]:
         return {
             "Oa": self.mo_coeff[0][:, :self.n_alpha],
             "Ob": self.mo_coeff[1][:, :self.n_beta],
@@ -189,7 +219,7 @@ def coefficients(self):
             "Vb": self.mo_coeff[1][:, self.n_beta:],
         }
 
-    def compute_mo_eri(self, blocks, spins):
+    def compute_mo_eri(self, blocks: Block4D, spins: Spin4D) -> Array4D:
         coeffs = tuple(self.coefficients[blocks[i] + spins[i]] for i in range(4))
         # TODO Pyscf uses HDF5 internal to do the AO2MO here we read it all
         #      into memory. This wastes memory and could be avoided if temporary
@@ -201,29 +231,30 @@ def compute_mo_eri(self, blocks, spins):
                                                     sizes[2], sizes[3])
 
 
-class PyScfHFProvider(HartreeFockProvider):
+class PyScfHFProvider(libadcc.HartreeFockProvider):
     """
         This implementation is only valid
         if no orbital reordering is required.
     """
-    def __init__(self, scfres):
+    def __init__(self, scfres: scf.hf.SCF):
         # Do not forget the next line,
         # otherwise weird errors result
         super().__init__()
-        self.scfres = scfres
+        self.scfres: scf.hf.SCF = scfres
         n_alpha, n_beta = scfres.mol.nelec
-        self.eri_builder = PyScfEriBuilder(self.scfres, self.n_orbs,
-                                           self.n_orbs_alpha, n_alpha,
-                                           n_beta, self.restricted)
-        self.operator_integral_provider = PyScfOperatorIntegralProvider(
-            self.scfres
+        self.eri_builder: PyScfEriBuilder = PyScfEriBuilder(
+            self.scfres, self.n_orbs, self.n_orbs_alpha,
+            n_alpha, n_beta, self.restricted
+        )
+        self.operator_integral_provider: PyScfOperatorIntegralProvider = (
+            PyScfOperatorIntegralProvider(self.scfres)
         )
         if not self.restricted:
             assert self.scfres.mo_coeff[0].shape[1] == \
                 self.scfres.mo_coeff[1].shape[1]
 
-        self.environment = None
-        self.environment_implementation = None
+        self.environment: Environment | None = None
+        self.environment_implementation: EnvironmentImplementation | None = None
         if hasattr(self.scfres, "with_solvent"):
             if hasattr(self.scfres.with_solvent, "cppe_state"):
                 self.environment = "pe"
@@ -232,22 +263,22 @@ def __init__(self, scfres):
                 self.environment = "pcm"
                 self.environment_implementation = "ddcosmo"
 
-    def pe_energy(self, dm, elec_only=True):
+    def pe_energy(self, dm: libadcc.Tensor, elec_only: bool = True) -> float:
         pe_state = self.scfres.with_solvent
         e_pe, _ = pe_state.kernel(dm.to_ndarray(), elec_only=elec_only)
-        return e_pe
+        return float(e_pe)
 
-    def pcm_energy(self, dm):
+    def pcm_energy(self, dm: libadcc.Tensor) -> float:
         # Since eps (dielectric constant) is the only solvent parameter
         # in pyscf and there is no solvent data available in the
         # program, the user needs to adjust scfres.with_solvent.eps
         # manually to the optical dielectric constant (if non
         # equilibrium solvation is desired).
         V_pcm = self.scfres.with_solvent._B_dot_x(dm.to_ndarray())
-        return np.einsum("uv,uv->", dm.to_ndarray(), V_pcm)
+        return float(np.einsum("uv,uv->", dm.to_ndarray(), V_pcm))
 
     @property
-    def excitation_energy_corrections(self):
+    def excitation_energy_corrections(self) -> dict[str, EnergyCorrection]:
         ret = []
         if self.environment == "pe":
             ptlr = EnergyCorrection(
@@ -269,17 +300,17 @@ def excitation_energy_corrections(self):
             ret.extend([ptlr])
         return {ec.name: ec for ec in ret}
 
-    def get_backend(self):
+    def get_backend(self) -> str:
         return "pyscf"
 
-    def get_conv_tol(self):
+    def get_conv_tol(self) -> float:
         if self.scfres.conv_tol_grad is None:
             conv_tol = self.scfres.conv_tol
         else:
             conv_tol = max(self.scfres.conv_tol, self.scfres.conv_tol_grad**2)
-        return conv_tol
+        return float(conv_tol)
 
-    def get_restricted(self):
+    def get_restricted(self) -> bool:
         if isinstance(self.scfres.mo_occ, list):
             restricted = len(self.scfres.mo_occ) < 2
         elif isinstance(self.scfres.mo_occ, np.ndarray):
@@ -289,27 +320,29 @@ def get_restricted(self):
                                    "not determine restricted / unrestricted.")
         return restricted
 
-    def get_energy_scf(self):
+    def get_energy_scf(self) -> float:
         return float(self.scfres.e_tot)
 
-    def get_nuclear_repulsion_energy(self):
+    def get_nuclear_repulsion_energy(self) -> float:
         return float(self.scfres.energy_nuc())
 
-    def get_spin_multiplicity(self):
+    def get_spin_multiplicity(self) -> int:
         # Note: In the pyscf world spin is 2S, so the multiplicity
         #       is spin + 1
         return int(self.scfres.mol.spin) + 1
 
-    def get_n_orbs_alpha(self):
+    def get_n_orbs_alpha(self) -> int:
         if self.restricted:
             return self.scfres.mo_coeff.shape[1]
         else:
             return self.scfres.mo_coeff[0].shape[1]
 
-    def get_n_bas(self):
+    def get_n_bas(self) -> int:
         return int(self.scfres.mol.nao_nr())
 
-    def get_nuclear_multipole(self, order, gauge_origin=(0, 0, 0)):
+    def get_nuclear_multipole(
+        self, order: int, gauge_origin: Coordinate = (0.0, 0.0, 0.0)
+    ) -> Array1D:
         charges = self.scfres.mol.atom_charges()
         if order == 0:
             # The function interface needs to be a np.array on return
@@ -326,10 +359,10 @@ def get_nuclear_multipole(self, order, gauge_origin=(0, 0, 0)):
         else:
             raise NotImplementedError("get_nuclear_multipole with order > 2")
 
-    def transform_gauge_origin_to_xyz(self, gauge_origin):
+    def transform_gauge_origin_to_xyz(self, gauge_origin: str) -> Coordinate:
         return _transform_gauge_origin_to_xyz(self.scfres, gauge_origin)
 
-    def fill_occupation_f(self, out):
+    def fill_occupation_f(self, out: Array1D) -> None:
         if self.restricted:
             out[:] = np.hstack((self.scfres.mo_occ / 2,
                                 self.scfres.mo_occ / 2))
@@ -337,7 +370,7 @@ def fill_occupation_f(self, out):
             out[:] = np.hstack((self.scfres.mo_occ[0],
                                 self.scfres.mo_occ[1]))
 
-    def fill_orbcoeff_fb(self, out):
+    def fill_orbcoeff_fb(self, out: Array2D) -> None:
         if self.restricted:
             mo_coeff = (self.scfres.mo_coeff,
                         self.scfres.mo_coeff)
@@ -347,7 +380,7 @@ def fill_orbcoeff_fb(self, out):
             np.hstack((mo_coeff[0], mo_coeff[1]))
         )
 
-    def fill_orben_f(self, out):
+    def fill_orben_f(self, out: Array1D) -> None:
         if self.restricted:
             out[:] = np.hstack((self.scfres.mo_energy,
                                 self.scfres.mo_energy))
@@ -355,25 +388,29 @@ def fill_orben_f(self, out):
             out[:] = np.hstack((self.scfres.mo_energy[0],
                                 self.scfres.mo_energy[1]))
 
-    def fill_fock_ff(self, slices, out):
+    def fill_fock_ff(self, slices: tuple[slice, slice], out: Array2D) -> None:
         diagonal = np.empty(self.n_orbs)
         self.fill_orben_f(diagonal)
         out[:] = np.diag(diagonal)[slices]
 
-    def fill_eri_ffff(self, slices, out):
+    def fill_eri_ffff(
+        self, slices: tuple[slice, slice, slice, slice], out: Array4D
+    ) -> None:
         self.eri_builder.fill_slice_symm(slices, out)
 
-    def fill_eri_phys_asym_ffff(self, slices, out):
+    def fill_eri_phys_asym_ffff(
+        self, slices: tuple[slice, slice, slice, slice], out: Array4D
+    ) -> None:
         raise NotImplementedError("fill_eri_phys_asym_ffff not implemented.")
 
-    def has_eri_phys_asym_ffff(self):
+    def has_eri_phys_asym_ffff(self) -> bool:
         return False
 
-    def flush_cache(self):
+    def flush_cache(self) -> None:
         self.eri_builder.flush_cache()
 
 
-def import_scf(scfres):
+def import_scf(scfres: scf.hf.SCF) -> PyScfHFProvider:
     # TODO The error messages here could be a bit more verbose
 
     if not isinstance(scfres, scf.hf.SCF):
@@ -391,8 +428,9 @@ def import_scf(scfres):
     return PyScfHFProvider(scfres)
 
 
-def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
-           conv_tol_grad=1e-9, max_iter=150, pe_options=None):
+def run_hf(xyz: str, basis: str, charge: int = 0, multiplicity: int = 1,
+           conv_tol: float = 1e-11, conv_tol_grad: float = 1e-9,
+           max_iter: int = 150, pe_options: dict | None = None) -> scf.hf.SCF:
     mol = gto.M(
         atom=xyz,
         basis=basis,
@@ -424,8 +462,9 @@ def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
     return mf
 
 
-def run_core_hole(xyz, basis, charge=0, multiplicity=1,
-                  conv_tol=1e-11, conv_tol_grad=1e-9, max_iter=150):
+def run_core_hole(xyz: str, basis: str, charge: int = 0, multiplicity: int = 1,
+                  conv_tol: float = 1e-11, conv_tol_grad: float = 1e-9,
+                  max_iter: int = 150) -> scf.hf.SCF:
     mol = gto.M(
         atom=xyz,
         basis=basis,
@@ -471,7 +510,9 @@ def run_core_hole(xyz, basis, charge=0, multiplicity=1,
     return mf_chole
 
 
-def _transform_gauge_origin_to_xyz(scfres, gauge_origin):
+def _transform_gauge_origin_to_xyz(
+    scfres: scf.hf.SCF, gauge_origin: Coordinate | str
+) -> Coordinate:
     """
     Determines the gauge origin. If the gauge origin is defined as a tuple
     the coordinates need to be given in atomic units!
@@ -492,6 +533,5 @@ def _transform_gauge_origin_to_xyz(scfres, gauge_origin):
         raise NotImplementedError("The gauge origin can be defined either by a "
                                   "keyword (origin, mass_center or charge_center) "
                                   "or by a tuple defining the Cartesian components "
-                                  "e.g. (x, y, z)."
-                                  )
+                                  "e.g. (x, y, z).")
     return gauge_origin
diff --git a/adcc/backends/veloxchem.py b/adcc/backends/veloxchem.py
index 4a99abd7..96b41195 100644
--- a/adcc/backends/veloxchem.py
+++ b/adcc/backends/veloxchem.py
@@ -96,9 +96,10 @@ def __init__(self, task, mol_orbs, n_orbs, n_orbs_alpha, n_alpha,
         super().__init__(n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted)
 
     def compute_mo_eri(self, blocks, spins):
-        eri = self.moints_drv.compute_in_memory(*self.compute_args,
-                                                moints_name="chem_" + blocks,
-                                                moints_spin=spins)
+        eri = self.moints_drv.compute_in_memory(
+            *self.compute_args, moints_name="chem_" + "".join(blocks),
+            moints_spin="".join(spins)
+        )
         return eri
 
 
@@ -171,6 +172,11 @@ def get_restricted(self):
     def get_energy_scf(self):
         return self.scfdrv.get_scf_energy()
 
+    def get_nuclear_repulsion_energy(self):
+        raise NotImplementedError(
+            "Nuclear repulsion energy not implemented for Veloxchem"
+        )
+
     def get_spin_multiplicity(self):
         return self.molecule.get_multiplicity()
 
diff --git a/conda/environment.yml b/conda/environment.yml
index 5de1d279..c375f2dd 100644
--- a/conda/environment.yml
+++ b/conda/environment.yml
@@ -8,7 +8,7 @@ dependencies:
   - gcc_linux-64
   - gxx_linux-64
   - libtensorlight >=3.0.1
-  - pybind11 >=2.6
+  - pybind11 >=3.0
   - pip
   - pkg-config
   # Run
diff --git a/libadcc.pyi b/libadcc.pyi
index 440820f8..54dc8c5a 100644
--- a/libadcc.pyi
+++ b/libadcc.pyi
@@ -1,6 +1,7 @@
 from __future__ import annotations
+import collections.abc
 import numpy
-import pybind11_stubgen.typing_ext
+import numpy.typing
 import typing
 
 __all__: list[str] = [
@@ -38,7 +39,10 @@ class AdcMemory:
     def __init__(self) -> None: ...
     def __repr__(self) -> str: ...
     def initialise(
-        self, pagefile_directory: str, max_block_size: int, allocator: str
+        self,
+        pagefile_directory: str,
+        max_block_size: typing.SupportsInt | typing.SupportsIndex = 16,
+        allocator: str = "standard",
     ) -> None: ...
     @property
     def allocator(self) -> str:
@@ -51,7 +55,9 @@ class AdcMemory:
         Get or set the batch size for contraction, i.e. the number of elements handled simultaneously in a tensor contraction.
         """
     @contraction_batch_size.setter
-    def contraction_batch_size(self, arg1: int) -> None: ...
+    def contraction_batch_size(
+        self, bsize: typing.SupportsInt | typing.SupportsIndex
+    ) -> None: ...
     @property
     def max_block_size(self) -> int:
         """
@@ -70,27 +76,45 @@ class HartreeFockProvider(HartreeFockSolution_i):
     In the remaining documentation we denote with `nf` the value returned by `get_n_orbs_alpha()` and with `nb` the value returned by `get_nbas()`.
     """
     def __init__(self) -> None: ...
-    def fill_eri_ffff(self, arg0: tuple, arg1: numpy.ndarray) -> None:
+    def fill_eri_ffff(
+        self,
+        slices: tuple[slice, slice, slice, slice],
+        out: numpy.ndarray[tuple[int, int, int, int], numpy.dtype[numpy.float64]],
+    ) -> None:
         """
-        Fill the passed numpy array `arg1` with a part of the electron-repulsion integral tensor in the molecular orbital basis. The indexing convention is the chemist's notation, i.e. the index tuple `(i,j,k,l)` refers to the integral :math:`(ij|kl)`. The block to store is specified by the provided tuple of ranges `arg0`, which gives the range of indices to place into the buffer along each of the axis. The index counting is done in spin orbitals, so the full range in each axis is `range(0, 2 * nf)`.
+        Fill the passed numpy array `out` with a part of the electron-repulsion integral tensor in the molecular orbital basis. The indexing convention is the chemist's notation, i.e. the index tuple `(i,j,k,l)` refers to the integral :math:`(ij|kl)`. The block to store is specified by the provided tuple of ranges `slices`, which gives the range of indices to place into the buffer along each of the axis. The index counting is done in spin orbitals, so the full range in each axis is `range(0, 2 * nf)`.
         """
-    def fill_eri_phys_asym_ffff(self, arg0: tuple, arg1: numpy.ndarray) -> None:
+    def fill_eri_phys_asym_ffff(
+        self,
+        slices: tuple[slice, slice, slice, slice],
+        out: numpy.ndarray[tuple[int, int, int, int], numpy.dtype[numpy.float64]],
+    ) -> None:
         """
-        Fill the passed numpy array `arg1` with a part of the **antisymmetrised** electron-repulsion integral tensor in the molecular orbital basis. The indexing convention is the physicist's notation, i.e. the index tuple `(i,j,k,l)` refers to the integral :math:`\\langle ij||kl \\rangle`. The block to store is specified by the provided tuple of ranges `arg0`, which gives the range of indices to place into the buffer along each of the axis. The index counting is done in spin orbitals, so the full range in each axis is `range(0, 2 * nf)`.
+        Fill the passed numpy array `out` with a part of the **antisymmetrised** electron-repulsion integral tensor in the molecular orbital basis. The indexing convention is the physicist's notation, i.e. the index tuple `(i,j,k,l)` refers to the integral :math:`\\langle ij||kl \\rangle`. The block to store is specified by the provided tuple of ranges `slices`, which gives the range of indices to place into the buffer along each of the axis. The index counting is done in spin orbitals, so the full range in each axis is `range(0, 2 * nf)`.
         """
-    def fill_fock_ff(self, arg0: tuple, arg1: numpy.ndarray) -> None:
+    def fill_fock_ff(
+        self,
+        slices: tuple[slice, slice],
+        out: numpy.ndarray[tuple[int, int], numpy.dtype[numpy.float64]],
+    ) -> None:
         """
-        Fill the passed numpy array `arg1` with a part of the Fock matrix in the molecular orbital basis. The block to store is specified by the provided tuple of ranges `arg0`, which gives the range of indices to place into the buffer along each of the axis. The index counting is done in spin orbitals, so the full range in each axis is `range(0, 2 * nf)`. The implementation should not assume that the alpha-beta and beta-alpha blocks are not accessed even though they are zero by spin symmetry.
+        Fill the passed numpy array `out` with a part of the Fock matrix in the molecular orbital basis. The block to store is specified by the provided tuple of ranges `slices`, which gives the range of indices to place into the buffer along each of the axis. The index counting is done in spin orbitals, so the full range in each axis is `range(0, 2 * nf)`. The implementation should not assume that the alpha-beta and beta-alpha blocks are not accessed even though they are zero by spin symmetry.
         """
-    def fill_occupation_f(self, arg0: numpy.ndarray) -> None:
+    def fill_occupation_f(
+        self, out: numpy.ndarray[tuple[int], numpy.dtype[numpy.float64]]
+    ) -> None:
         """
         Fill the passed numpy array of size `(2 * nf, )` with the occupation number for each SCF orbital.
         """
-    def fill_orbcoeff_fb(self, arg0: numpy.ndarray) -> None:
+    def fill_orbcoeff_fb(
+        self, out: numpy.ndarray[tuple[int, int], numpy.dtype[numpy.float64]]
+    ) -> None:
         """
         Fill the passed numpy array of size `(2 * nf, nb)` with the SCF orbital coefficients, i.e. the uniform transform from the one-particle basis to the molecular orbitals.
         """
-    def fill_orben_f(self, arg0: numpy.ndarray) -> None:
+    def fill_orben_f(
+        self, out: numpy.ndarray[tuple[int], numpy.dtype[numpy.float64]]
+    ) -> None:
         """
         Fill the passed numpy array of size `(2 * nf, )` with the SCF orbital energies.
         """
@@ -120,8 +144,8 @@ class HartreeFockProvider(HartreeFockSolution_i):
         Returns the number of HF *spin* orbitals of alpha spin. It is assumed the same number of beta spin orbitals are used. This value is abbreviated by `nf` in the documentation.
         """
     def get_nuclear_multipole(
-        self, arg0: int, arg1: tuple
-    ) -> numpy.ndarray[numpy.float64]:
+        self, order: int, gauge_origin: tuple[float, float, float] = (0.0, 0.0, 0.0)
+    ) -> numpy.ndarray[tuple[int], numpy.dtype[numpy.float64]]:
         """
         Returns the nuclear multipole of the requested order. For `0` returns the total nuclear charge as an array of size 1, for `1` returns the nuclear dipole moment as an array of size 3.
         """
@@ -141,7 +165,9 @@ class HartreeFockProvider(HartreeFockSolution_i):
         """
         Returns whether `fill_eri_phys_asym_ffff` function is implemented and should be used(*True*) or whether antisymmetrisation should be done inside adcc starting from the `fill_eri_ffff` function (*False*)
         """
-    def transform_gauge_origin_to_xyz(self, arg0: str) -> tuple:
+    def transform_gauge_origin_to_xyz(
+        self, gauge_origin: str
+    ) -> tuple[float, float, float]:
         """
         Transforms a string specifying the gauge origin to a tuple containing the x, y, z Cartesian components.
         """
@@ -157,7 +183,7 @@ class HartreeFockSolution_i:
     @property
     def energy_scf(self) -> float: ...
     @property
-    def fock_ff(self) -> numpy.ndarray[numpy.float64]: ...
+    def fock_ff(self) -> numpy.ndarray[tuple[int, int], numpy.dtype[numpy.float64]]: ...
     @property
     def n_alpha(self) -> int: ...
     @property
@@ -173,11 +199,13 @@ class HartreeFockSolution_i:
     @property
     def nuclear_repulsion_energy(self) -> float: ...
     @property
-    def occupation_f(self) -> numpy.ndarray[numpy.float64]: ...
+    def occupation_f(self) -> numpy.ndarray[tuple[int], numpy.dtype[numpy.float64]]: ...
     @property
-    def orbcoeff_fb(self) -> numpy.ndarray[numpy.float64]: ...
+    def orbcoeff_fb(
+        self,
+    ) -> numpy.ndarray[tuple[int, int], numpy.dtype[numpy.float64]]: ...
     @property
-    def orben_f(self) -> numpy.ndarray[numpy.float64]: ...
+    def orben_f(self) -> numpy.ndarray[tuple[int], numpy.dtype[numpy.float64]]: ...
     @property
     def restricted(self) -> bool: ...
     @property
@@ -188,20 +216,26 @@ class MoIndexTranslation:
     Helper object to extract information from indices into orbitals subspaces and to map them between different indexing conventions (full MO space, MO subspaces, indexing convention in the HF Provider / SCF host program, ... Python binding to :cpp:class:`libadcc::MoIndexTranslation`.
     """
     @typing.overload
-    def __init__(self, arg0: MoSpaces, arg1: str) -> None:
+    def __init__(self, mospaces: MoSpaces, space: str) -> None:
         """
         Construct a MoIndexTranslation class from an MoSpaces object and the identifier for the space (e.g. o1o1, v1o1, o3v2o1v1, ...)
         """
     @typing.overload
-    def __init__(self, arg0: MoSpaces, arg1: list[str]) -> None:
+    def __init__(
+        self, mospaces: MoSpaces, subspaces: collections.abc.Sequence[str]
+    ) -> None:
         """
         Construct a MoIndexTranslation class from an MoSpaces object and the list of identifiers for the space (e.g. ["o1", "o1"] ...)
         """
-    def block_index_of(self, arg0: tuple) -> tuple:
+    def block_index_of(
+        self, tpl: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> tuple[int, ...]:
         """
         Get the block index of an index, i.e. get the index which points to the block of the tensor in which the element with the passed index is contained in.
         """
-    def block_index_spatial_of(self, arg0: tuple) -> tuple:
+    def block_index_spatial_of(
+        self, tpl: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> tuple[int, ...]:
         """
         Get the spatial block index of an index
 
@@ -212,36 +246,65 @@ class MoIndexTranslation:
         map to the same value upon a call of this function.
         """
     @typing.overload
-    def combine(self, arg0: tuple, arg1: tuple) -> tuple:
+    def combine(
+        self,
+        bidx: tuple[typing.SupportsInt | typing.SupportsIndex, ...],
+        ibidx: tuple[typing.SupportsInt | typing.SupportsIndex, ...],
+    ) -> tuple[int, ...]:
         """
         Combine a block index and an in-block index into the appropriate index. Effectively undoes the effect of 'split'.
         """
     @typing.overload
-    def combine(self, arg0: str, arg1: tuple, arg2: tuple) -> tuple:
+    def combine(
+        self,
+        spin_block: str,
+        bidx: tuple[typing.SupportsInt | typing.SupportsIndex, ...],
+        ibidx: tuple[typing.SupportsInt | typing.SupportsIndex, ...],
+    ) -> tuple[int, ...]:
         """
         Combine a spin block (given as a string of 'a's or 'b's), a spatial-only block index and an in-block index into the appropriate index. Essentially undoes the effect of 'spin_of', 'block_index_spatial_of' and 'inblock_index_of'.
         """
-    def full_index_of(self, arg0: tuple) -> tuple:
+    def full_index_of(
+        self, tpl: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> tuple[int, ...]:
         """
         Map an index given in the space, which was passed upon construction, to the corresponding index in the full MO index range (the ffff space).
         """
-    def hf_provider_index_of(self, arg0: tuple) -> tuple:
+    def hf_provider_index_of(
+        self, index: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> tuple[int, ...]:
         """
         Map an index (given in the space passed upon construction) to the indexing convention of the host program provided to adcc as the HF provider.
         """
-    def inblock_index_of(self, arg0: tuple) -> tuple:
+    def inblock_index_of(
+        self, tpl: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> tuple[int, ...]:
         """
         Get the in-block index, i.e. the index within the tensor block.
         """
-    def map_range_to_hf_provider(self, arg0: tuple) -> list:
+    def map_range_to_hf_provider(
+        self,
+        ranges: tuple[
+            tuple[
+                typing.SupportsInt | typing.SupportsIndex,
+                typing.SupportsInt | typing.SupportsIndex,
+            ],
+            ...,
+        ],
+    ) -> list[dict[str, tuple[tuple[int, int], ...]]]:
         """
         Map a range of indices to host program indices, i.e. the indexing convention
         used in the HfProvider, which provides the SCF data to adcc.
 
         Since the mapping between subspace and host program indices might not be contiguous,
-        a list of pairs of ranges is returned. In each pair, the first entry represents a
-        range of indices (indexed in the MO subspace) and the second entry represents
-        the equivalent range of indices in the Hartree-Fock provider these are mapped to.
+        a list of range mappings is returned. Each mapping is a dict with the keys
+        'from' and 'to': 'from' holds a range of indices (indexed in the MO subspace)
+        and 'to' the equivalent range of indices in the Hartree-Fock provider these are
+        mapped to. Both are given as a tuple of half-open [start, end) index pairs, one
+        pair for each dimension, i.e. in the same format as the `ranges` argument.
+
+        For example, for a two-dimensional space:
+          [{'from': ((0, 1), (0, 1)), 'to': ((0, 1), (2, 3))}, ...]
 
           ranges    Tuple of pairs of indices: One index pair for each dimension. Each
                     pair describes the range of indices along one axis, which should be
@@ -249,15 +312,21 @@ class MoIndexTranslation:
                     be thought of as a half-open interval [start, end), where start and
                     end are the indexed passed as a pair to the function.
         """
-    def spin_of(self, arg0: tuple) -> str:
+    def spin_of(
+        self, tpl: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> str:
         """
         Get the spin block of each of the index components as a string.
         """
-    def split(self, arg0: tuple) -> tuple:
+    def split(
+        self, tpl: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> tuple[tuple[int, ...], tuple[int, ...]]:
         """
         Split an index into block index and in-block index
         """
-    def split_spin(self, arg0: tuple) -> tuple:
+    def split_spin(
+        self, tpl: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> tuple[str, tuple[int, ...], tuple[int, ...]]:
         """
         Split an index into a spin block descriptor, a spatial block index and an in-block index.
         """
@@ -272,7 +341,7 @@ class MoIndexTranslation:
         Return the number of dimensions.
         """
     @property
-    def shape(self) -> tuple:
+    def shape(self) -> tuple[int, ...]:
         """
         Return the length along each dimension.
         """
@@ -282,7 +351,10 @@ class MoIndexTranslation:
         Return the space supplied on initialisation.
         """
     @property
-    def subspaces(self) -> list[str]: ...
+    def subspaces(self) -> list[str]:
+        """
+        Return the space supplied on initialisation split into the subspace along each dimension, e.g. ["o1", "v1"] for the space "o1v1".
+        """
 
 class MoSpaces:
     """
@@ -290,17 +362,23 @@ class MoSpaces:
     """
     def __init__(
         self,
-        arg0: HartreeFockSolution_i,
-        arg1: AdcMemory,
-        arg2: list[int],
-        arg3: list[int],
-        arg4: list[int],
+        hf: HartreeFockSolution_i,
+        adcmem: AdcMemory,
+        core_orbitals: collections.abc.Sequence[
+            typing.SupportsInt | typing.SupportsIndex
+        ],
+        frozen_core_orbitals: collections.abc.Sequence[
+            typing.SupportsInt | typing.SupportsIndex
+        ],
+        frozen_virtuals: collections.abc.Sequence[
+            typing.SupportsInt | typing.SupportsIndex
+        ],
     ) -> None:
         """
         Construct an MoSpaces object from a HartreeFockSolution_i, a pointer to
         an AdcMemory object.
 
-        adcmem_ptr        ADC memory keep-alive object to be used in all Tensors
+        adcmem            ADC memory keep-alive object to be used in all Tensors
                           constructed using this MoSpaces object.
         core_orbitals     List of orbitals indices (in the full fock space, original
                           ordering of the hf object), which defines the orbitals to
@@ -317,15 +395,15 @@ class MoSpaces:
                           in the ADC calculation. The same number of alpha and beta
                           orbitals has to be selected.
         """
-    def n_orbs(self, arg0: str) -> int:
+    def n_orbs(self, space: str) -> int:
         """
         The number of orbitals in a particular orbital subspace
         """
-    def n_orbs_alpha(self, arg0: str) -> int:
+    def n_orbs_alpha(self, space: str) -> int:
         """
         The number of alpha orbitals in a particular orbital subspace
         """
-    def n_orbs_beta(self, arg0: str) -> int:
+    def n_orbs_beta(self, space: str) -> int:
         """
         The number of beta orbitals in a particular orbital subspace
         """
@@ -397,14 +475,19 @@ class ReferenceState:
     """
     Class representing information about the reference state for adcc. Python binding to:cpp:class:`libadcc::ReferenceState`.
     """
-    def __init__(self, arg0: HartreeFockSolution_i, arg1: MoSpaces, arg2: bool) -> None:
+    def __init__(
+        self,
+        hfsoln: HartreeFockSolution_i,
+        mo: MoSpaces,
+        symmetry_check_on_import: bool,
+    ) -> None:
         """
         Setup a ReferenceStateject using an MoSpaces object.
 
-        hfsoln_ptr        Pointer to the Interface to the host program,
+        hfsoln            Pointer to the Interface to the host program,
                           providing the HartreeFockSolution data, which
                           will be provided by this object.
-        mo_ptr            MoSpaces object containing info about the MoSpace setup
+        mo                MoSpaces object containing info about the MoSpace setup
                           and the point group symmetry.
         symmetry_check_on_import
                           Should symmetry of the imported objects be checked
@@ -414,7 +497,7 @@ class ReferenceState:
                           from the host programs. Do not enable this unless you know
                           that you really want to.
         """
-    def eri(self, arg0: str) -> Tensor:
+    def eri(self, space: str) -> Tensor:
         """
         Return the ERI (electron-repulsion integrals) tensor block corresponding to the provided space.
         """
@@ -422,32 +505,36 @@ class ReferenceState:
         """
         Tell the contained HartreeFockSolution_i object (which was passed upon construction), that a larger amount of import operations is done and that the next request for further imports will most likely take some time, such that intermediate caches can now be flushed to save some memory or other resources.
         """
-    def fock(self, arg0: str) -> Tensor:
+    def fock(self, space: str) -> Tensor:
         """
         Return the Fock matrix block corresponding to the provided space.
         """
-    def gauge_origin_to_xyz(self, arg0: str) -> tuple: ...
+    def gauge_origin_to_xyz(self, gauge_origin: str) -> tuple[float, float, float]: ...
     def import_all(self) -> None:
         """
         Normally the class only imports the Fock matrix blocks and electron-repulsion integrals of a particular space combination when this is requested by a call to above fock() or eri() functions. This function call, however, instructs the class to immediately import *all* such blocks. Typically you do not want to do this.
         """
     def nuclear_quadrupole(
         self,
-        arg0: typing.Annotated[list[float], pybind11_stubgen.typing_ext.FixedSize(3)],
-    ) -> numpy.ndarray[numpy.float64]: ...
-    def orbital_coefficients(self, arg0: str) -> Tensor:
+        gauge_origin: tuple[
+            typing.SupportsFloat | typing.SupportsIndex,
+            typing.SupportsFloat | typing.SupportsIndex,
+            typing.SupportsFloat | typing.SupportsIndex,
+        ],
+    ) -> numpy.ndarray[tuple[int], numpy.dtype[numpy.float64]]: ...
+    def orbital_coefficients(self, space: str) -> Tensor:
         """
         Return the molecular orbital coefficients corresponding to the provided space (alpha and beta coefficients are returned)
         """
-    def orbital_coefficients_alpha(self, arg0: str) -> Tensor:
+    def orbital_coefficients_alpha(self, space: str) -> Tensor:
         """
         Return the alpha molecular orbital coefficients corresponding to the provided space
         """
-    def orbital_coefficients_beta(self, arg0: str) -> Tensor:
+    def orbital_coefficients_beta(self, space: str) -> Tensor:
         """
         Return the beta molecular orbital coefficients corresponding to the provided space
         """
-    def orbital_energies(self, arg0: str) -> Tensor:
+    def orbital_energies(self, space: str) -> Tensor:
         """
         Return the orbital energies corresponding to the provided space
         """
@@ -464,7 +551,7 @@ class ReferenceState:
         Setting this property allows to drop ERI tensor blocks if they are no longer needed to save memory.
         """
     @cached_eri_blocks.setter
-    def cached_eri_blocks(self, arg1: list[str]) -> None: ...
+    def cached_eri_blocks(self, newlist: collections.abc.Sequence[str]) -> None: ...
     @property
     def cached_fock_blocks(self) -> list[str]:
         """
@@ -473,7 +560,7 @@ class ReferenceState:
         Setting this property allows to drop fock matrix blocks if they are no longer needed to save memory.
         """
     @cached_fock_blocks.setter
-    def cached_fock_blocks(self, arg1: list[str]) -> None: ...
+    def cached_fock_blocks(self, newlist: collections.abc.Sequence[str]) -> None: ...
     @property
     def conv_tol(self) -> float:
         """
@@ -525,7 +612,9 @@ class ReferenceState:
         Number of beta orbitals
         """
     @property
-    def nuclear_dipole(self) -> numpy.ndarray[numpy.float64]: ...
+    def nuclear_dipole(
+        self,
+    ) -> numpy.ndarray[tuple[int], numpy.dtype[numpy.float64]]: ...
     @property
     def nuclear_repulsion_energy(self) -> float:
         """
@@ -554,13 +643,22 @@ class Symmetry:
     Container for Tensor symmetry information
     """
     @typing.overload
-    def __init__(self, arg0: MoSpaces, arg1: str) -> None:
+    def __init__(self, mospaces: MoSpaces, space: str) -> None:
         """
         Construct a Symmetry class from an MoSpaces object and the identifier for the space (e.g. o1o1, v1o1, o3v2o1v1, ...). Python binding to :cpp:class:`libadcc::Symmetry`.
         """
     @typing.overload
     def __init__(
-        self, arg0: MoSpaces, arg1: str, arg2: dict[str, tuple[int, int]]
+        self,
+        mospaces: MoSpaces,
+        space: str,
+        extra_axes_orbs: collections.abc.Mapping[
+            str,
+            tuple[
+                typing.SupportsInt | typing.SupportsIndex,
+                typing.SupportsInt | typing.SupportsIndex,
+            ],
+        ],
     ) -> None:
         """
         Construct a Symmetry class from an MoSpaces object, a space string and a map to supply the number of orbitals for some additional axes.
@@ -587,7 +685,7 @@ class Symmetry:
         The list of irreducible representations, for which the tensor shall be non-zero. If this is *not* set, i.e. an empty list, all irreps will be allowed.
         """
     @irreps_allowed.setter
-    def irreps_allowed(self, arg1: list[str]) -> None: ...
+    def irreps_allowed(self, irreps: collections.abc.Sequence[str]) -> None: ...
     @property
     def mospaces(self) -> MoSpaces:
         """
@@ -612,9 +710,9 @@ class Symmetry:
         is only rudimentary at the moment.
         """
     @permutations.setter
-    def permutations(self, arg1: list[str]) -> None: ...
+    def permutations(self, permutations: collections.abc.Sequence[str]) -> None: ...
     @property
-    def shape(self) -> tuple:
+    def shape(self) -> tuple[int, ...]:
         """
         Return the shape of tensors constructed from this symmetry.
         """
@@ -631,7 +729,12 @@ class Symmetry:
         with a factor of -1.0 between them.
         """
     @spin_block_maps.setter
-    def spin_block_maps(self, arg1: list[tuple[str, str, float]]) -> None: ...
+    def spin_block_maps(
+        self,
+        spin_maps: collections.abc.Sequence[
+            tuple[str, str, typing.SupportsFloat | typing.SupportsIndex]
+        ],
+    ) -> None: ...
     @property
     def spin_blocks_forbidden(self) -> list[str]:
         """
@@ -639,72 +742,100 @@ class Symmetry:
         Blocks are given as a string in the letters 'a' and 'b', e.g. ["aaba", "abba"]
         """
     @spin_blocks_forbidden.setter
-    def spin_blocks_forbidden(self, arg1: list[str]) -> None: ...
+    def spin_blocks_forbidden(
+        self, forbidden: collections.abc.Sequence[str]
+    ) -> None: ...
 
 class Tensor:
     """
     Class representing the Tensor objects used for computations in adcc
     """
-
-    flags: list[str]
     @typing.overload
-    def __add__(self, arg0: float) -> Tensor: ...
+    def __add__(
+        self, number: typing.SupportsFloat | typing.SupportsIndex
+    ) -> Tensor: ...
     @typing.overload
-    def __add__(self, arg0: Tensor) -> Tensor: ...
-    def __getitem__(self, arg0: tuple) -> float:
+    def __add__(self, other: Tensor) -> Tensor: ...
+    def __getitem__(
+        self, idcs: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> float:
         """
         Get a tensor element or a slice of tensor elements.
         """
-    def __iadd__(self, arg0: Tensor) -> Tensor: ...
-    def __imul__(self, arg0: float) -> Tensor: ...
-    def __init__(self, arg0: Symmetry) -> None:
+    def __iadd__(self, other: Tensor) -> Tensor: ...
+    def __imul__(
+        self, number: typing.SupportsFloat | typing.SupportsIndex
+    ) -> Tensor: ...
+    def __init__(self, symmetry: Symmetry) -> None:
         """
         Construct a Tensor object using a Symmetry object describing its symmetry properties.
         The returned object is not guaranteed to contain initialised memory. Python binding to :cpp:class:`libadcc::Tensor`
         """
-    def __isub__(self, arg0: Tensor) -> Tensor: ...
-    def __itruediv__(self, arg0: float) -> Tensor: ...
+    def __isub__(self, other: Tensor) -> Tensor: ...
+    def __itruediv__(
+        self, number: typing.SupportsFloat | typing.SupportsIndex
+    ) -> Tensor: ...
     def __len__(self) -> int: ...
-    def __matmul__(self, arg0: Tensor) -> Tensor: ...
+    def __matmul__(self, other: Tensor) -> Tensor: ...
     @typing.overload
-    def __mul__(self, arg0: float) -> Tensor: ...
+    def __mul__(
+        self, number: typing.SupportsFloat | typing.SupportsIndex
+    ) -> Tensor: ...
     @typing.overload
-    def __mul__(self, arg0: Tensor) -> Tensor:
+    def __mul__(self, other: Tensor) -> Tensor:
         """
         Multiply two tensors elementwise.
         """
     def __neg__(self) -> Tensor: ...
     def __pos__(self) -> Tensor: ...
-    def __radd__(self, arg0: float) -> Tensor: ...
-    def __repr__(self) -> typing.Any: ...
-    def __rmul__(self, arg0: float) -> Tensor: ...
-    def __rsub__(self, arg0: float) -> Tensor: ...
-    def __setitem__(self, arg0: tuple, arg1: float) -> float:
+    def __radd__(
+        self, number: typing.SupportsFloat | typing.SupportsIndex
+    ) -> Tensor: ...
+    def __repr__(self) -> str: ...
+    def __rmul__(
+        self, number: typing.SupportsFloat | typing.SupportsIndex
+    ) -> Tensor: ...
+    def __rsub__(
+        self, number: typing.SupportsFloat | typing.SupportsIndex
+    ) -> Tensor: ...
+    def __setitem__(
+        self,
+        idcs: tuple[typing.SupportsInt | typing.SupportsIndex, ...],
+        value: typing.SupportsFloat | typing.SupportsIndex,
+    ) -> None:
         """
         Set a tensor element or a slice of tensor elements. The operation will adhere symmetry, i.e. alter all elements equivalent by symmetry at once.
         """
-    def __str__(self) -> typing.Any: ...
+    def __str__(self) -> str: ...
     @typing.overload
-    def __sub__(self, arg0: float) -> Tensor: ...
+    def __sub__(
+        self, number: typing.SupportsFloat | typing.SupportsIndex
+    ) -> Tensor: ...
     @typing.overload
-    def __sub__(self, arg0: Tensor) -> Tensor: ...
+    def __sub__(self, other: Tensor) -> Tensor: ...
     @typing.overload
-    def __truediv__(self, arg0: float) -> Tensor: ...
+    def __truediv__(
+        self, number: typing.SupportsFloat | typing.SupportsIndex
+    ) -> Tensor: ...
     @typing.overload
-    def __truediv__(self, arg0: Tensor) -> Tensor:
+    def __truediv__(self, other: Tensor) -> Tensor:
         """
         Divide two tensors elementwise.
         """
     @typing.overload
-    def antisymmetrise(self, arg0: list) -> Tensor: ...
+    def antisymmetrise(
+        self,
+        permutations: collections.abc.Iterable[int]
+        | collections.abc.Iterable[collections.abc.Iterable[int]],
+    ) -> Tensor: ...
     @typing.overload
-    def antisymmetrise(self, *args) -> Tensor: ...
+    def antisymmetrise(self, *args: int) -> Tensor: ...
     def copy(self) -> Tensor:
         """
         Returns a deep copy of the tensor.
         """
     @typing.overload
-    def describe_expression(self, arg0: str) -> str:
+    def describe_expression(self, stage: str) -> str:
         """
         Return a string providing a hopefully descriptive representation of the tensor expression stored inside the object.
         """
@@ -714,46 +845,62 @@ class Tensor:
         """
         Return a string providing a hopefully descriptive representation of the symmetry information stored inside the tensor.
         """
-    def diagonal(self, *args) -> Tensor: ...
+    def diagonal(self, *args: int) -> Tensor: ...
     @typing.overload
-    def dot(self, arg0: Tensor) -> float: ...
+    def dot(self, other: Tensor) -> float: ...
     @typing.overload
-    def dot(self, arg0: list) -> numpy.ndarray[numpy.float64]: ...
+    def dot(
+        self, tensors: list[Tensor]
+    ) -> numpy.ndarray[tuple[int], numpy.dtype[numpy.float64]]: ...
     def empty_like(self) -> Tensor: ...
     def evaluate(self) -> Tensor:
         """
         Ensure the tensor to be fully evaluated and resilient in memory. Usually happens automatically when needed. Might be useful for fine-tuning, however.
         """
-    def is_allowed(self, arg0: tuple) -> bool:
+    def is_allowed(
+        self, idcs: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> bool:
         """
         Is a particular index allowed by symmetry
         """
     def nosym_like(self) -> Tensor: ...
     def ones_like(self) -> Tensor: ...
-    def select_n_absmax(self, arg0: int) -> list:
+    def select_n_absmax(
+        self, n: typing.SupportsInt | typing.SupportsIndex
+    ) -> list[tuple[list[int], float]]:
         """
         Select the n absolute maximal elements.
         """
-    def select_n_absmin(self, arg0: int) -> list:
+    def select_n_absmin(
+        self, n: typing.SupportsInt | typing.SupportsIndex
+    ) -> list[tuple[list[int], float]]:
         """
         Select the n absolute minimal elements.
         """
-    def select_n_max(self, arg0: int) -> list:
+    def select_n_max(
+        self, n: typing.SupportsInt | typing.SupportsIndex
+    ) -> list[tuple[list[int], float]]:
         """
         Select the n maximal elements.
         """
-    def select_n_min(self, arg0: int) -> list:
+    def select_n_min(
+        self, n: typing.SupportsInt | typing.SupportsIndex
+    ) -> list[tuple[list[int], float]]:
         """
         Select the n minimal elements.
         """
     @typing.overload
-    def set_from_ndarray(self, arg0: numpy.ndarray) -> Tensor:
+    def set_from_ndarray(
+        self, in_array: typing.Annotated[numpy.typing.ArrayLike, numpy.float64]
+    ) -> Tensor:
         """
-        Set all tensor elements from a standard np::ndarray by making a copy. Provide an optional tolerance argument to increase the tolerance for the check for symmetry consistency.
+        Set all tensor elements from a standard np::ndarray by making a copy.
         """
     @typing.overload
     def set_from_ndarray(
-        self, arg0: numpy.ndarray[numpy.float64], arg1: float
+        self,
+        in_array: typing.Annotated[numpy.typing.ArrayLike, numpy.float64],
+        symmetry_tolerance: typing.SupportsFloat | typing.SupportsIndex,
     ) -> Tensor:
         """
         Set all tensor elements from a standard np::ndarray by making a copy. Provide an optional tolerance argument to increase the tolerance for the check for symmetry consistency.
@@ -762,7 +909,9 @@ class Tensor:
         """
         Set the tensor as immutable, allowing some optimisations to be performed.
         """
-    def set_mask(self, arg0: str, arg1: float) -> None:
+    def set_mask(
+        self, mask: str, value: typing.SupportsFloat | typing.SupportsIndex
+    ) -> None:
         """
         Set all elements corresponding to an index mask, which is given by a string eg. 'iijkli' sets elements T_{iijkli}
         """
@@ -771,21 +920,31 @@ class Tensor:
         Set all tensor elements to random data, adhering to the internal symmetry.
         """
     @typing.overload
-    def symmetrise(self, arg0: list) -> Tensor: ...
+    def symmetrise(
+        self,
+        permutations: collections.abc.Iterable[int]
+        | collections.abc.Iterable[collections.abc.Iterable[int]],
+    ) -> Tensor: ...
     @typing.overload
-    def symmetrise(self, *args) -> Tensor: ...
-    def to_ndarray(self) -> numpy.ndarray[numpy.float64]:
+    def symmetrise(self, *args: int) -> Tensor: ...
+    def to_ndarray(self) -> numpy.typing.NDArray[numpy.float64]:
         """
         Export the tensor data to a standard np::ndarray by making a copy.
         """
     @typing.overload
     def transpose(self) -> Tensor: ...
     @typing.overload
-    def transpose(self, arg0: tuple) -> Tensor: ...
+    def transpose(
+        self, axes: tuple[typing.SupportsInt | typing.SupportsIndex, ...]
+    ) -> Tensor: ...
     def zeros_like(self) -> Tensor: ...
     @property
     def T(self) -> Tensor: ...
     @property
+    def flags(self) -> list[str]: ...
+    @flags.setter
+    def flags(self, new_flags: collections.abc.Sequence[str]) -> None: ...
+    @property
     def mutable(self) -> bool: ...
     @property
     def ndim(self) -> int: ...
@@ -795,7 +954,7 @@ class Tensor:
         Does the tensor need evaluation or is it fully evaluated and resilient in memory.
         """
     @property
-    def shape(self) -> tuple: ...
+    def shape(self) -> tuple[int, ...]: ...
     @property
     def size(self) -> int: ...
     @property
@@ -803,28 +962,30 @@ class Tensor:
     @property
     def subspaces(self) -> list[str]: ...
 
-def amplitude_vector_enforce_spin_kind(arg0: Tensor, arg1: str, arg2: str) -> None:
+def amplitude_vector_enforce_spin_kind(
+    doubles_tensor: Tensor, block: str, spin_kind: str
+) -> None:
     """
     Apply the spin symmetrisation required to make the doubles and higher parts of an amplitude vector consist of components for a particular spin kind only.
     """
 
 def direct_sum(a: Tensor, b: Tensor) -> Tensor: ...
-def evaluate(arg0: Tensor) -> Tensor: ...
+def evaluate(tensor: Tensor) -> Tensor: ...
 def fill_pp_doubles_guesses(
-    guesses_d: list[Tensor],
+    guesses_d: collections.abc.Sequence[Tensor],
     mospaces: MoSpaces,
-    df02: Tensor,
-    df13: Tensor,
-    spin_change_twice: int,
-    degeneracy_tolerance: float,
+    df1: Tensor,
+    df2: Tensor,
+    spin_change_twice: typing.SupportsInt | typing.SupportsIndex,
+    degeneracy_tolerance: typing.SupportsFloat | typing.SupportsIndex,
 ) -> int:
     """
-    Fill the passed vector of doubles blocks with doubles guesses using the delta-Fock matrices df02 and df13, which are the two delta-Fock matrices involved in the doubles block.
+    Fill the passed vector of doubles blocks with doubles guesses using the delta-Fock matrices df1 and df2, which are the two delta-Fock matrices involved in the doubles block.
 
     guesses_d    Vectors of guesses, all elements are assumed to be initialised to zero and the symmetry is assumed to be properly set up.
     mospaces     Mospaces object
-    df02         Delta-Fock between spaces 0 and 2 of the ADC matrix
-    df13         Delta-Fock between spaces 1 and 3 of the ADC matrix
+    df1         Delta-Fock between spaces 0 and 2 of the ADC matrix
+    df2         Delta-Fock between spaces 1 and 3 of the ADC matrix
     spin_change_twice   Twice the value of the spin change to enforce in an excitation.
     degeneracy_tolerance  Tolerance for two entries of the diagonal to be considered degenerate, i.e. identical.
     Returns     The number of guess vectors which have been properly initialised (the others are invalid and should be discarded).
@@ -841,9 +1002,10 @@ def get_n_threads_total() -> int:
     """
 
 def linear_combination_strict(
-    coefficients: numpy.ndarray[numpy.float64], tensors: list
+    coefficients: typing.Annotated[numpy.typing.ArrayLike, numpy.float64],
+    tensors: list[Tensor],
 ) -> Tensor: ...
-def make_symmetry_eri(arg0: MoSpaces, arg1: str) -> Symmetry:
+def make_symmetry_eri(mospaces: MoSpaces, space: str) -> Symmetry:
     """
     Return the Symmetry object like it would be set up for the passed subspace
     of the electron-repulsion tensor.
@@ -852,17 +1014,23 @@ def make_symmetry_eri(arg0: MoSpaces, arg1: str) -> Symmetry:
       space       Space string (e.g. o1v1o1v1)
     """
 
-def make_symmetry_operator(arg0: MoSpaces, arg1: str, arg2: str, arg3: str) -> Symmetry:
+def make_symmetry_operator(
+    mospaces: MoSpaces,
+    space: str,
+    operator_symmetry: str,
+    cartesian_transformation: str,
+) -> Symmetry:
     """
     Return the Symmetry object for an orbital subspace block of a one-particle operator
 
       mospaces    MoSpaces object
       space       Space string (e.g. o1v1)
-      symmetry    Describes the symmetry of the tensor (only in effect if both
+      operator_symmetry
+                  Describes the symmetry of the tensor (only in effect if both
                   subspaces of the space string are identical).
                   Valid are "nosymmetry", "hermitian" and "antihermitian".
       cartesian_transformation
-                  The cartesian function according to which the operator transforms.
+                  The cartesian function according to which the operator               transforms.
 
     Valid cartesian_transformation values include:
          "1"                   Totally symmetric (default)
@@ -872,7 +1040,11 @@ def make_symmetry_operator(arg0: MoSpaces, arg1: str, arg2: str, arg3: str) -> S
     """
 
 def make_symmetry_operator_basis(
-    arg0: MoSpaces, arg1: int, arg2: str, arg3: int, arg4: str
+    mospaces: MoSpaces,
+    n_bas: typing.SupportsInt | typing.SupportsIndex,
+    operator_symmetry: str,
+    n_particle_op: typing.SupportsInt | typing.SupportsIndex,
+    blocks: str,
 ) -> Symmetry:
     """
     Return the symmetry object for an operator in the AO basis. The object will
@@ -882,11 +1054,11 @@ def make_symmetry_operator_basis(
     where M is an n_bas x n_bas block and is indentical in upper-left
     and lower-right.
 
-    mospaces_ptr      MoSpaces pointer
+    mospaces          MoSpaces pointer
     n_bas             Number of AO basis functions
     operator_symmetry Is the tensor symmetric (hermitian/antihermitian, only
                       in effect if both space axes identical).
-                      Nosymmetry disables a setup of permutational symmetry.
+                      Nosymmetry disables a setup of 'bra-ket' symmetry.
     n_particle_op     NParticleOperator
     blocks            Which blocks of the operator to return. Valid values
                       are 'ab' to return a tensor for both alpha and beta
@@ -895,7 +1067,10 @@ def make_symmetry_operator_basis(
     """
 
 def make_symmetry_orbital_coefficients(
-    arg0: MoSpaces, arg1: str, arg2: int, arg3: str
+    mospaces: MoSpaces,
+    space: str,
+    n_bas: typing.SupportsInt | typing.SupportsIndex,
+    blocks: str = "ab",
 ) -> Symmetry:
     """
     Return the Symmetry object like it would be set up for the passed subspace
@@ -907,7 +1082,7 @@ def make_symmetry_orbital_coefficients(
       blocks      Spin blocks to include. Valid are "ab", "a" and "b".
     """
 
-def make_symmetry_orbital_energies(arg0: MoSpaces, arg1: str) -> Symmetry:
+def make_symmetry_orbital_energies(mospaces: MoSpaces, space: str) -> Symmetry:
     """
     Return the Symmetry object like it would be set up for the passed subspace
     of the orbital energies tensor.
@@ -916,7 +1091,7 @@ def make_symmetry_orbital_energies(arg0: MoSpaces, arg1: str) -> Symmetry:
       space       space string (e.g. o1)
     """
 
-def make_symmetry_triples(arg0: MoSpaces, arg1: str) -> Symmetry:
+def make_symmetry_triples(mospaces: MoSpaces, space: str) -> Symmetry:
     """
     Return the Symmetry object like it would be set up for the passed subspace
     of a triples amplitude tensor.
@@ -925,22 +1100,26 @@ def make_symmetry_triples(arg0: MoSpaces, arg1: str) -> Symmetry:
       space       Space string (e.g. o1o1o1v1v1v1)
     """
 
-def set_n_threads(arg0: int) -> None:
+def set_n_threads(n_threads: typing.SupportsInt | typing.SupportsIndex) -> None:
     """
     Set the number of running worker threads used by adcc
     """
 
-def set_n_threads_total(arg0: int) -> None:
+def set_n_threads_total(n_total: typing.SupportsInt | typing.SupportsIndex) -> None:
     """
     Set the total number of threads (running and sleeping) used by adcc. This will disappear in the future. Do not rely on it.
     """
 
 @typing.overload
-def tensordot(a: Tensor, b: Tensor, axes: typing.Iterable) -> typing.Any: ...
+def tensordot(
+    a: Tensor, b: Tensor, axes: collections.abc.Iterable[collections.abc.Iterable[int]]
+) -> Tensor | float: ...
 @typing.overload
-def tensordot(a: Tensor, b: Tensor, axes: int) -> typing.Any: ...
+def tensordot(
+    a: Tensor, b: Tensor, axes: typing.SupportsInt | typing.SupportsIndex
+) -> Tensor | float: ...
 @typing.overload
-def tensordot(a: Tensor, b: Tensor) -> typing.Any: ...
+def tensordot(a: Tensor, b: Tensor) -> Tensor | float: ...
 @typing.overload
 def trace(subscripts: str, tensor: Tensor) -> float: ...
 @typing.overload
diff --git a/libadcc_src/HartreeFockSolution_i.hh b/libadcc_src/HartreeFockSolution_i.hh
index 914523ee..f6a0292e 100644
--- a/libadcc_src/HartreeFockSolution_i.hh
+++ b/libadcc_src/HartreeFockSolution_i.hh
@@ -46,12 +46,13 @@ class HartreeFockSolution_i {
   ///@{
   /** Fill a buffer with nuclear multipole data for the nuclear multipole of
    *  given order. */
-  virtual void nuclear_multipole(size_t order, std::array<scalar_type, 3> gauge_origin,
-                                 scalar_type* buffer, size_t size) const = 0;
+  virtual void nuclear_multipole(
+        size_t order, std::tuple<scalar_type, scalar_type, scalar_type> gauge_origin,
+        scalar_type* buffer, size_t size) const = 0;
   //@}
 
   /** Determine the gauge origin. */
-  virtual const std::array<scalar_type, 3> gauge_origin_to_xyz(
+  virtual const std::tuple<scalar_type, scalar_type, scalar_type> gauge_origin_to_xyz(
         std::string gauge_origin) const = 0;
 
   /** \name Sizes of the data */
diff --git a/libadcc_src/ReferenceState.cc b/libadcc_src/ReferenceState.cc
index 7e32326a..4d680112 100644
--- a/libadcc_src/ReferenceState.cc
+++ b/libadcc_src/ReferenceState.cc
@@ -19,7 +19,6 @@
 
 #include "ReferenceState.hh"
 #include "MoIndexTranslation.hh"
-#include "TensorImpl.hh"
 #include "exceptions.hh"
 #include "import_eri.hh"
 #include "make_symmetry.hh"
@@ -351,16 +350,16 @@ std::string ReferenceState::irreducible_representation() const {
 }
 
 std::vector<scalar_type> ReferenceState::nuclear_multipole(
-      size_t order, std::array<scalar_type, 3> gauge_origin) const {
+      size_t order,
+      std::tuple<scalar_type, scalar_type, scalar_type> gauge_origin) const {
   std::vector<scalar_type> ret((order + 2) * (order + 1) / 2);
   m_hfsoln_ptr->nuclear_multipole(order, gauge_origin, ret.data(), ret.size());
   return ret;
 }
 
-const std::array<scalar_type, 3> ReferenceState::gauge_origin_to_xyz(
-      std::string gauge_origin) const {
-  const std::array<scalar_type, 3> ret = m_hfsoln_ptr->gauge_origin_to_xyz(gauge_origin);
-  return ret;
+const std::tuple<scalar_type, scalar_type, scalar_type>
+ReferenceState::gauge_origin_to_xyz(std::string gauge_origin) const {
+  return m_hfsoln_ptr->gauge_origin_to_xyz(gauge_origin);
 }
 
 //
diff --git a/libadcc_src/ReferenceState.hh b/libadcc_src/ReferenceState.hh
index 72f0115d..621cac05 100644
--- a/libadcc_src/ReferenceState.hh
+++ b/libadcc_src/ReferenceState.hh
@@ -18,7 +18,6 @@
 //
 
 #pragma once
-#include "AdcMemory.hh"
 #include "HartreeFockSolution_i.hh"
 #include "MoSpaces.hh"
 #include "Tensor.hh"
@@ -92,12 +91,13 @@ class ReferenceState {
 
   /** Return the nuclear contribution to the cartesian multipole moment
    *  (in standard ordering, i.e. xx, xy, xz, yy, yz, zz) of the given order. */
-  std::vector<scalar_type> nuclear_multipole(size_t order,
-                                             std::array<scalar_type, 3> gauge_origin = {
-                                                   0, 0, 0}) const;
+  std::vector<scalar_type> nuclear_multipole(
+        size_t order,
+        std::tuple<scalar_type, scalar_type, scalar_type> gauge_origin = {0, 0, 0}) const;
 
   /** Determine the gauge origin for nuclear multipoles. */
-  const std::array<scalar_type, 3> gauge_origin_to_xyz(std::string gauge_origin) const;
+  const std::tuple<scalar_type, scalar_type, scalar_type> gauge_origin_to_xyz(
+        std::string gauge_origin) const;
 
   /** Return the SCF convergence tolerance */
   double conv_tol() const { return m_hfsoln_ptr->conv_tol(); }
diff --git a/libadcc_src/make_symmetry.hh b/libadcc_src/make_symmetry.hh
index 84dee8b6..99269691 100644
--- a/libadcc_src/make_symmetry.hh
+++ b/libadcc_src/make_symmetry.hh
@@ -95,7 +95,7 @@ std::shared_ptr<Symmetry> make_symmetry_operator(
  * \param n_bas               Number of AO basis functions
  * \param operator_symmetry   Is the tensor symmetric (hermitian/antihermitian, only
  *                            in effect if both space axes identical).
- *                            Nosymmetry disables a setup of permutational symmetry.
+ *                            Nosymmetry disables a setup of 'bra-ket' symmetry.
  * \param n_particle_op       NParticle Operator
  * \param blocks              Which blocks of the operator to return. Valid values
  *                            are "ab" to return a tensor for both alpha and beta
diff --git a/libadcc_src/pyiface/export_AdcMemory.cc b/libadcc_src/pyiface/export_AdcMemory.cc
index ba74582e..c330500c 100644
--- a/libadcc_src/pyiface/export_AdcMemory.cc
+++ b/libadcc_src/pyiface/export_AdcMemory.cc
@@ -24,7 +24,6 @@
 
 namespace libadcc {
 
-using namespace pybind11::literals;
 namespace py = pybind11;
 
 static std::string AdcMemory___repr__(const AdcMemory& self) {
@@ -40,11 +39,11 @@ static std::string AdcMemory___repr__(const AdcMemory& self) {
 }
 
 void export_AdcMemory(py::module& m) {
-  py::class_<AdcMemory, std::shared_ptr<AdcMemory>>(
+  py::class_<AdcMemory, std::shared_ptr<AdcMemory>> adc_memory(
         m, "AdcMemory",
         "Class controlling the memory allocations for adcc ADC calculations. Python "
-        "binding to :cpp:class:`libadcc::AdcMemory`.")
-        .def(py::init<>())
+        "binding to :cpp:class:`libadcc::AdcMemory`.");
+  adc_memory.def(py::init<>())
         .def_property_readonly("allocator", &AdcMemory::allocator,
                                "Return the allocator to which the class is initialised.")
         .def_property_readonly("pagefile_directory", &AdcMemory::pagefile_directory,
@@ -54,11 +53,12 @@ void export_AdcMemory(py::module& m) {
               "max_block_size", &AdcMemory::max_block_size,
               "Return the maximal block size a tenor may have along each axis.")
         .def_property("contraction_batch_size", &AdcMemory::contraction_batch_size,
-                      &AdcMemory::set_contraction_batch_size,
+                      py::cpp_function(&AdcMemory::set_contraction_batch_size,
+                                       py::is_method(adc_memory), py::arg("bsize")),
                       "Get or set the batch size for contraction, i.e. the number of "
                       "elements handled simultaneously in a tensor contraction.")
-        .def("initialise", &AdcMemory::initialise, "pagefile_directory"_a,
-             "max_block_size"_a, "allocator"_a)
+        .def("initialise", &AdcMemory::initialise, py::arg("pagefile_directory"),
+             py::arg("max_block_size") = 16, py::arg("allocator") = "standard")
         .def("__repr__", &AdcMemory___repr__)
         //
         ;
diff --git a/libadcc_src/pyiface/export_HartreeFockProvider.cc b/libadcc_src/pyiface/export_HartreeFockProvider.cc
index 2acc5fc7..b9187552 100644
--- a/libadcc_src/pyiface/export_HartreeFockProvider.cc
+++ b/libadcc_src/pyiface/export_HartreeFockProvider.cc
@@ -20,6 +20,7 @@
 #include "../HartreeFockSolution_i.hh"
 #include "../exceptions.hh"
 #include "hartree_fock_solution_hack.hh"
+#include "ndarray.hh"
 #include "util.hh"
 #include <pybind11/pybind11.h>
 
@@ -50,9 +51,10 @@ class HartreeFockProvider : public HartreeFockSolution_i {
   //
   // Translate C++-like interface to python-like interface
   //
-  void nuclear_multipole(size_t order, std::array<scalar_type, 3> gauge_origin,
+  void nuclear_multipole(size_t order,
+                         std::tuple<scalar_type, scalar_type, scalar_type> gauge_origin,
                          scalar_type* buffer, size_t size) const override {
-    py::array_t<scalar_type> ret = get_nuclear_multipole(order, py::cast(gauge_origin));
+    py::array_t<scalar_type> ret = get_nuclear_multipole(order, gauge_origin);
     if (static_cast<ssize_t>(size) != ret.size()) {
       throw dimension_mismatch("Array size (==" + std::to_string(ret.size()) +
                                ") does not agree with buffer size (" +
@@ -61,10 +63,9 @@ class HartreeFockProvider : public HartreeFockSolution_i {
     std::copy(ret.data(), ret.data() + size, buffer);
   }
 
-  const std::array<scalar_type, 3> gauge_origin_to_xyz(
+  const std::tuple<scalar_type, scalar_type, scalar_type> gauge_origin_to_xyz(
         std::string gauge_origin) const override {
-    return py::cast<std::array<scalar_type, 3>>(
-          transform_gauge_origin_to_xyz(py::cast(gauge_origin)));
+    return transform_gauge_origin_to_xyz(py::cast(gauge_origin));
   }
 
   void occupation_f(scalar_type* buffer, size_t size) const override {
@@ -133,9 +134,9 @@ class HartreeFockProvider : public HartreeFockSolution_i {
           buffer, d1_length * d2_length * sizeof(scalar_type));
     std::vector<ssize_t> strides{static_cast<ssize_t>(sizeof(scalar_type) * d1_stride),
                                  static_cast<ssize_t>(sizeof(scalar_type) * d2_stride)};
-    py::tuple slices = py::make_tuple(
+    const std::tuple<py::slice, py::slice> slices{
           py::slice(static_cast<ssize_t>(d1_start), static_cast<ssize_t>(d1_end), 1),
-          py::slice(static_cast<ssize_t>(d2_start), static_cast<ssize_t>(d2_end), 1));
+          py::slice(static_cast<ssize_t>(d2_start), static_cast<ssize_t>(d2_end), 1)};
     fill_fock_ff(slices, py::array({d1_length, d2_length}, strides, buffer, memview));
   }
 
@@ -186,11 +187,11 @@ class HartreeFockProvider : public HartreeFockSolution_i {
                                  static_cast<ssize_t>(sizeof(scalar_type) * d2_stride),
                                  static_cast<ssize_t>(sizeof(scalar_type) * d3_stride),
                                  static_cast<ssize_t>(sizeof(scalar_type) * d4_stride)};
-    py::tuple slices = py::make_tuple(
+    const std::tuple<py::slice, py::slice, py::slice, py::slice> slices{
           py::slice(static_cast<ssize_t>(d1_start), static_cast<ssize_t>(d1_end), 1),
           py::slice(static_cast<ssize_t>(d2_start), static_cast<ssize_t>(d2_end), 1),
           py::slice(static_cast<ssize_t>(d3_start), static_cast<ssize_t>(d3_end), 1),
-          py::slice(static_cast<ssize_t>(d4_start), static_cast<ssize_t>(d4_end), 1));
+          py::slice(static_cast<ssize_t>(d4_start), static_cast<ssize_t>(d4_end), 1)};
     fill_eri_ffff(slices, py::array({d1_length, d2_length, d3_length, d4_length}, strides,
                                     buffer, memview));
   }
@@ -243,11 +244,11 @@ class HartreeFockProvider : public HartreeFockSolution_i {
                                  static_cast<ssize_t>(sizeof(scalar_type) * d2_stride),
                                  static_cast<ssize_t>(sizeof(scalar_type) * d3_stride),
                                  static_cast<ssize_t>(sizeof(scalar_type) * d4_stride)};
-    py::tuple slices = py::make_tuple(
+    std::tuple<py::slice, py::slice, py::slice, py::slice> slices{
           py::slice(static_cast<ssize_t>(d1_start), static_cast<ssize_t>(d1_end), 1),
           py::slice(static_cast<ssize_t>(d2_start), static_cast<ssize_t>(d2_end), 1),
           py::slice(static_cast<ssize_t>(d3_start), static_cast<ssize_t>(d3_end), 1),
-          py::slice(static_cast<ssize_t>(d4_start), static_cast<ssize_t>(d4_end), 1));
+          py::slice(static_cast<ssize_t>(d4_start), static_cast<ssize_t>(d4_end), 1)};
     fill_eri_phys_asym_ffff(
           slices, py::array({d1_length, d2_length, d3_length, d4_length}, strides, buffer,
                             memview));
@@ -258,22 +259,29 @@ class HartreeFockProvider : public HartreeFockSolution_i {
   //
   virtual size_t get_n_orbs_alpha() const = 0;
   virtual size_t get_n_bas() const        = 0;
-  virtual py::array_t<scalar_type> get_nuclear_multipole(
-        size_t order, py::tuple gauge_origin) const                                 = 0;
-  virtual const py::tuple transform_gauge_origin_to_xyz(py::str gauge_origin) const = 0;
-  virtual real_type get_conv_tol() const                                            = 0;
-  virtual bool get_restricted() const                                               = 0;
-  virtual size_t get_spin_multiplicity() const                                      = 0;
-  virtual real_type get_energy_scf() const                                          = 0;
-  virtual real_type get_nuclear_repulsion_energy() const                            = 0;
-  virtual std::string get_backend() const                                           = 0;
-
-  virtual void fill_occupation_f(py::array out) const                         = 0;
-  virtual void fill_orben_f(py::array out) const                              = 0;
-  virtual void fill_orbcoeff_fb(py::array out) const                          = 0;
-  virtual void fill_fock_ff(py::tuple, py::array out) const                   = 0;
-  virtual void fill_eri_ffff(py::tuple slices, py::array out) const           = 0;
-  virtual void fill_eri_phys_asym_ffff(py::tuple slices, py::array out) const = 0;
+  virtual NDArray<scalar_type, 1> get_nuclear_multipole(
+        py::int_ order,
+        std::tuple<py::float_, py::float_, py::float_> gauge_origin) const = 0;
+  virtual const std::tuple<scalar_type, scalar_type, scalar_type>
+  transform_gauge_origin_to_xyz(py::str gauge_origin) const = 0;
+  virtual real_type get_conv_tol() const                    = 0;
+  virtual bool get_restricted() const                       = 0;
+  virtual size_t get_spin_multiplicity() const              = 0;
+  virtual real_type get_energy_scf() const                  = 0;
+  virtual real_type get_nuclear_repulsion_energy() const    = 0;
+  virtual std::string get_backend() const                   = 0;
+
+  virtual void fill_occupation_f(NDArray<scalar_type, 1> out) const = 0;
+  virtual void fill_orben_f(NDArray<scalar_type, 1> out) const      = 0;
+  virtual void fill_orbcoeff_fb(NDArray<scalar_type, 2> out) const  = 0;
+  virtual void fill_fock_ff(std::tuple<py::slice, py::slice>,
+                            NDArray<scalar_type, 2> out) const      = 0;
+  virtual void fill_eri_ffff(
+        std::tuple<py::slice, py::slice, py::slice, py::slice> slices,
+        NDArray<scalar_type, 4> out) const = 0;
+  virtual void fill_eri_phys_asym_ffff(
+        std::tuple<py::slice, py::slice, py::slice, py::slice> slices,
+        NDArray<scalar_type, 4> out) const = 0;
 };
 
 /** This implements the trampoline for C++ to call the python functions
@@ -289,14 +297,17 @@ class PyHartreeFockProvider : public HartreeFockProvider {
   size_t get_n_bas() const override {
     PYBIND11_OVERLOAD_PURE(size_t, HartreeFockProvider, get_n_bas, );
   }
-  py::array_t<scalar_type> get_nuclear_multipole(size_t order,
-                                                 py::tuple gauge_origin) const override {
-    PYBIND11_OVERLOAD_PURE(py::array_t<scalar_type>, HartreeFockProvider,
+  NDArray<scalar_type, 1> get_nuclear_multipole(
+        py::int_ order,
+        std::tuple<py::float_, py::float_, py::float_> gauge_origin) const override {
+    PYBIND11_OVERLOAD_PURE(PYBIND11_TYPE(NDArray<scalar_type, 1>), HartreeFockProvider,
                            get_nuclear_multipole, order, gauge_origin);
   }
-  const py::tuple transform_gauge_origin_to_xyz(py::str gauge_origin) const override {
-    PYBIND11_OVERLOAD_PURE(py::tuple, HartreeFockProvider, transform_gauge_origin_to_xyz,
-                           gauge_origin);
+  const std::tuple<scalar_type, scalar_type, scalar_type> transform_gauge_origin_to_xyz(
+        py::str gauge_origin) const override {
+    PYBIND11_OVERLOAD_PURE(
+          PYBIND11_TYPE(std::tuple<scalar_type, scalar_type, scalar_type>),
+          HartreeFockProvider, transform_gauge_origin_to_xyz, gauge_origin);
   }
   real_type get_conv_tol() const override {
     PYBIND11_OVERLOAD_PURE(real_type, HartreeFockProvider, get_conv_tol, );
@@ -314,22 +325,26 @@ class PyHartreeFockProvider : public HartreeFockProvider {
     PYBIND11_OVERLOAD_PURE(real_type, HartreeFockProvider,
                            get_nuclear_repulsion_energy, );
   }
-  void fill_occupation_f(py::array out) const override {
+  void fill_occupation_f(NDArray<scalar_type, 1> out) const override {
     PYBIND11_OVERLOAD_PURE(void, HartreeFockProvider, fill_occupation_f, out);
   }
-  void fill_orben_f(py::array out) const override {
+  void fill_orben_f(NDArray<scalar_type, 1> out) const override {
     PYBIND11_OVERLOAD_PURE(void, HartreeFockProvider, fill_orben_f, out);
   }
-  void fill_orbcoeff_fb(py::array out) const override {
+  void fill_orbcoeff_fb(NDArray<scalar_type, 2> out) const override {
     PYBIND11_OVERLOAD_PURE(void, HartreeFockProvider, fill_orbcoeff_fb, out);
   }
-  void fill_fock_ff(py::tuple slices, py::array out) const override {
+  void fill_fock_ff(std::tuple<py::slice, py::slice> slices,
+                    NDArray<scalar_type, 2> out) const override {
     PYBIND11_OVERLOAD_PURE(void, HartreeFockProvider, fill_fock_ff, slices, out);
   }
-  void fill_eri_ffff(py::tuple slices, py::array out) const override {
+  void fill_eri_ffff(std::tuple<py::slice, py::slice, py::slice, py::slice> slices,
+                     NDArray<scalar_type, 4> out) const override {
     PYBIND11_OVERLOAD_PURE(void, HartreeFockProvider, fill_eri_ffff, slices, out);
   }
-  void fill_eri_phys_asym_ffff(py::tuple slices, py::array out) const override {
+  void fill_eri_phys_asym_ffff(
+        std::tuple<py::slice, py::slice, py::slice, py::slice> slices,
+        NDArray<scalar_type, 4> out) const override {
     PYBIND11_OVERLOAD_PURE(void, HartreeFockProvider, fill_eri_phys_asym_ffff, slices,
                            out);
   }
@@ -344,7 +359,7 @@ class PyHartreeFockProvider : public HartreeFockProvider {
   }
 };
 
-static py::array_t<scalar_type> HartreeFockSolution_i_occupation_f(
+static NDArray<scalar_type, 1> HartreeFockSolution_i_occupation_f(
       const HartreeFockSolution_i& self) {
   py::array_t<scalar_type> ret(self.n_orbs());
   self.occupation_f(ret.mutable_data(), self.n_orbs());
@@ -369,21 +384,21 @@ static size_t count_electrons(const HartreeFockSolution_i& self, bool count_beta
   return ret;
 }
 
-static py::array_t<scalar_type> HartreeFockSolution_i_orben_f(
+static NDArray<scalar_type, 1> HartreeFockSolution_i_orben_f(
       const HartreeFockSolution_i& self) {
   py::array_t<scalar_type> ret(self.n_orbs());
   self.orben_f(ret.mutable_data(), self.n_orbs());
   return ret;
 }
 
-static py::array_t<scalar_type> HartreeFockSolution_i_orbcoeff_fb(
+static NDArray<scalar_type, 2> HartreeFockSolution_i_orbcoeff_fb(
       const HartreeFockSolution_i& self) {
   py::array_t<scalar_type> ret({self.n_orbs(), self.n_bas()});
   self.orbcoeff_fb(ret.mutable_data(), self.n_orbs() * self.n_bas());
   return ret;
 }
 
-static py::array_t<scalar_type> HartreeFockSolution_i_fock_ff(
+static NDArray<scalar_type, 2> HartreeFockSolution_i_fock_ff(
       const HartreeFockSolution_i& self) {
   py::array_t<scalar_type> ret({self.n_orbs(), self.n_orbs()});
   self.fock_ff(0, self.n_orbs(), 0, self.n_orbs(),
@@ -461,50 +476,54 @@ void export_HartreeFockProvider(py::module& m) {
              "Returns the number of *spatial* one-electron basis functions. This value "
              "is abbreviated by `nb` in the documentation.")
         .def("get_nuclear_multipole", &HartreeFockProvider::get_nuclear_multipole,
+             py::arg("order"), py::arg("gauge_origin") = py::make_tuple(0.0, 0.0, 0.0),
              "Returns the nuclear multipole of the requested order. For `0` returns the "
              "total nuclear charge as an array of size 1, for `1` returns the nuclear "
              "dipole moment as an array of size 3.")
         //
         .def("transform_gauge_origin_to_xyz",
-             &HartreeFockProvider::transform_gauge_origin_to_xyz,
+             &HartreeFockProvider::transform_gauge_origin_to_xyz, py::arg("gauge_origin"),
              "Transforms a string specifying the gauge origin to a tuple containing "
              "the x, y, z Cartesian components.")
         //
-        .def("fill_occupation_f", &HartreeFockProvider::fill_occupation_f,
+        .def("fill_occupation_f", &HartreeFockProvider::fill_occupation_f, py::arg("out"),
              "Fill the passed numpy array of size `(2 * nf, )` with the occupation "
              "number for each SCF orbital.")
-        .def("fill_orben_f", &HartreeFockProvider::fill_orben_f,
+        .def("fill_orben_f", &HartreeFockProvider::fill_orben_f, py::arg("out"),
              "Fill the passed numpy array of size `(2 * nf, )` with the SCF orbital "
              "energies.")
-        .def("fill_orbcoeff_fb", &HartreeFockProvider::fill_orbcoeff_fb,
+        .def("fill_orbcoeff_fb", &HartreeFockProvider::fill_orbcoeff_fb, py::arg("out"),
              "Fill the passed numpy array of size `(2 * nf, nb)` with the SCF orbital "
              "coefficients, i.e. the uniform transform from the one-particle basis to "
              "the molecular orbitals.")
-        .def("fill_fock_ff", &HartreeFockProvider::fill_fock_ff,
-             "Fill the passed numpy array `arg1` with a part of the Fock matrix in the "
+        .def("fill_fock_ff", &HartreeFockProvider::fill_fock_ff, py::arg("slices"),
+             py::arg("out"),
+             "Fill the passed numpy array `out` with a part of the Fock matrix in the "
              "molecular orbital basis. The block to store is specified by the provided "
-             "tuple of ranges `arg0`, which gives the range of indices to place into the "
-             "buffer along each of the axis. The index counting is done in spin "
+             "tuple of ranges `slices`, which gives the range of indices to place into "
+             "the buffer along each of the axis. The index counting is done in spin "
              "orbitals, so the full range in each axis is `range(0, 2 * nf)`. The "
              "implementation should not assume that the alpha-beta and beta-alpha blocks "
              "are not accessed even though they are zero by spin symmetry.")
-        .def("fill_eri_ffff", &HartreeFockProvider::fill_eri_ffff,
-             "Fill the passed numpy array `arg1` with a part of the electron-repulsion "
+        .def("fill_eri_ffff", &HartreeFockProvider::fill_eri_ffff, py::arg("slices"),
+             py::arg("out"),
+             "Fill the passed numpy array `out` with a part of the electron-repulsion "
              "integral tensor in the molecular orbital basis. "
              "The indexing convention is the chemist's notation, i.e. the index tuple "
              "`(i,j,k,l)` refers to the integral :math:`(ij|kl)`. "
              "The block to store is specified by the provided "
-             "tuple of ranges `arg0`, which gives the range of indices to place into the "
-             "buffer along each of the axis. The index counting is done in spin "
+             "tuple of ranges `slices`, which gives the range of indices to place into "
+             "the buffer along each of the axis. The index counting is done in spin "
              "orbitals, so the full range in each axis is `range(0, 2 * nf)`.")
         .def("fill_eri_phys_asym_ffff", &HartreeFockProvider::fill_eri_phys_asym_ffff,
-             "Fill the passed numpy array `arg1` with a part of the **antisymmetrised** "
+             py::arg("slices"), py::arg("out"),
+             "Fill the passed numpy array `out` with a part of the **antisymmetrised** "
              "electron-repulsion integral tensor in the molecular orbital basis. "
              "The indexing convention is the physicist's notation, i.e. the index tuple "
              "`(i,j,k,l)` refers to the integral :math:`\\langle ij||kl \\rangle`. "
              "The block to store is specified by the provided "
-             "tuple of ranges `arg0`, which gives the range of indices to place into the "
-             "buffer along each of the axis. The index counting is done in spin "
+             "tuple of ranges `slices`, which gives the range of indices to place into "
+             "the buffer along each of the axis. The index counting is done in spin "
              "orbitals, so the full range in each axis is `range(0, 2 * nf)`.")
         .def("has_eri_phys_asym_ffff", &HartreeFockProvider::has_eri_phys_asym_ffff,
              "Returns whether `fill_eri_phys_asym_ffff` function is implemented and "
diff --git a/libadcc_src/pyiface/export_MoIndexTranslation.cc b/libadcc_src/pyiface/export_MoIndexTranslation.cc
index 8339c92b..7c68b3b6 100644
--- a/libadcc_src/pyiface/export_MoIndexTranslation.cc
+++ b/libadcc_src/pyiface/export_MoIndexTranslation.cc
@@ -21,11 +21,17 @@
 #include "util.hh"
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
+#include <pybind11/typing.h>
 
 namespace libadcc {
 
 namespace py = pybind11;
 
+// Type definition used throughout the interface
+using IdxTuple   = py::typing::Tuple<py::ssize_t, py::ellipsis>;
+using RangePair  = py::typing::Tuple<py::ssize_t, py::ssize_t>;
+using RangeTuple = py::typing::Tuple<RangePair, py::ellipsis>;
+
 static std::vector<size_t> parse_tuple(size_t ndim, const py::tuple& tuple) {
   if (tuple.size() != ndim) {
     throw py::value_error(
@@ -40,7 +46,7 @@ static std::vector<size_t> parse_tuple(size_t ndim, const py::tuple& tuple) {
   return ret;
 }
 
-static py::tuple convert_range_to_tuples(const SimpleRange& range) {
+static RangeTuple convert_range_to_tuples(const SimpleRange& range) {
   py::tuple ret(range.size());
   for (size_t i = 0; i < range.size(); ++i) {
     ret[i] = py::make_tuple(range[i].first, range[i].second);
@@ -48,8 +54,9 @@ static py::tuple convert_range_to_tuples(const SimpleRange& range) {
   return ret;
 }
 
-static py::list MoIndexTranslation_map_range_to_hf_provider(
-      std::shared_ptr<MoIndexTranslation> self, py::tuple ranges) {
+static py::typing::List<py::typing::Dict<py::str, RangeTuple>>
+MoIndexTranslation_map_range_to_hf_provider(std::shared_ptr<MoIndexTranslation> self,
+                                            RangeTuple ranges) {
   if (ranges.size() != self->ndim()) {
     throw py::value_error("Number of elements passed in the index range tuple (== " +
                           std::to_string(ranges.size()) + ") and dimensionality (== " +
@@ -97,13 +104,18 @@ void export_MoIndexTranslation(py::module& m) {
         "subspaces, indexing convention in the HF Provider / SCF host program, ... "
         "Python binding to :cpp:class:`libadcc::MoIndexTranslation`.")
         .def(py::init<std::shared_ptr<const MoSpaces>, const std::string&>(),
+             py::arg("mospaces"), py::arg("space"),
              "Construct a MoIndexTranslation class from an MoSpaces object and the "
              "identifier for "
              "the space (e.g. o1o1, v1o1, o3v2o1v1, ...)")
         .def(py::init<std::shared_ptr<const MoSpaces>, const std::vector<std::string>&>(),
+             py::arg("mospaces"), py::arg("subspaces"),
              "Construct a MoIndexTranslation class from an MoSpaces object and the "
              "list of identifiers for the space (e.g. [\"o1\", \"o1\"] ...)")
-        .def_property_readonly("subspaces", &MoIndexTranslation::subspaces)
+        .def_property_readonly("subspaces", &MoIndexTranslation::subspaces,
+                               "Return the space supplied on initialisation split into "
+                               "the subspace along each dimension, e.g. [\"o1\", \"v1\"] "
+                               "for the space \"o1v1\".")
         .def_property_readonly("mospaces", &MoIndexTranslation::mospaces_ptr,
                                "Return the MoSpaces object supplied on initialisation")
         .def_property_readonly("space", &MoIndexTranslation::space,
@@ -119,25 +131,28 @@ void export_MoIndexTranslation(py::module& m) {
         //
         .def(
               "full_index_of",
-              [](std::shared_ptr<MoIndexTranslation> self, py::tuple tpl) {
+              [](std::shared_ptr<MoIndexTranslation> self, IdxTuple tpl) {
                 return shape_tuple(self->full_index_of(parse_tuple(self->ndim(), tpl)));
               },
+              py::arg("tpl"),
               "Map an index given in the space, which was passed upon construction, to "
               "the corresponding index in the full MO index range (the ffff space).")
         .def(
               "block_index_of",
-              [](std::shared_ptr<MoIndexTranslation> self, py::tuple tpl) {
+              [](std::shared_ptr<MoIndexTranslation> self, IdxTuple tpl) {
                 return shape_tuple(self->block_index_of(parse_tuple(self->ndim(), tpl)));
               },
+              py::arg("tpl"),
               "Get the block index of an index, i.e. get the index which points to the "
               "block of the tensor in which the element with the passed index is "
               "contained in.")
         .def(
               "block_index_spatial_of",
-              [](std::shared_ptr<MoIndexTranslation> self, py::tuple tpl) {
+              [](std::shared_ptr<MoIndexTranslation> self, IdxTuple tpl) {
                 return shape_tuple(
                       self->block_index_spatial_of(parse_tuple(self->ndim(), tpl)));
               },
+              py::arg("tpl"),
               "Get the spatial block index of an index\n"
               "\n"
               "The spatial block index is the result of block_index_of modulo the spin "
@@ -150,79 +165,95 @@ void export_MoIndexTranslation(py::module& m) {
               "map to the same value upon a call of this function.")
         .def(
               "inblock_index_of",
-              [](std::shared_ptr<MoIndexTranslation> self, py::tuple tpl) {
+              [](std::shared_ptr<MoIndexTranslation> self, IdxTuple tpl) {
                 return shape_tuple(
                       self->inblock_index_of(parse_tuple(self->ndim(), tpl)));
               },
+              py::arg("tpl"),
               "Get the in-block index, i.e. the index within the tensor block.")
         .def(
               "spin_of",
-              [](std::shared_ptr<MoIndexTranslation> self, py::tuple tpl) {
+              [](std::shared_ptr<MoIndexTranslation> self, IdxTuple tpl) {
                 return self->spin_of(parse_tuple(self->ndim(), tpl));
               },
+              py::arg("tpl"),
               "Get the spin block of each of the index components as a string.")
         .def(
               "split",
-              [](std::shared_ptr<MoIndexTranslation> self, py::tuple tpl) {
+              [](std::shared_ptr<MoIndexTranslation> self,
+                 IdxTuple tpl) -> py::typing::Tuple<IdxTuple, IdxTuple> {
                 auto splitted = self->split(parse_tuple(self->ndim(), tpl));
                 return py::make_tuple(shape_tuple(splitted.first),
                                       shape_tuple(splitted.second));
               },
-              "Split an index into block index and in-block index")
+              py::arg("tpl"), "Split an index into block index and in-block index")
         .def(
               "split_spin",
-              [](std::shared_ptr<MoIndexTranslation> self, py::tuple tpl) {
+              [](std::shared_ptr<MoIndexTranslation> self,
+                 IdxTuple tpl) -> py::typing::Tuple<py::str, IdxTuple, IdxTuple> {
                 auto splitted = self->split_spin(parse_tuple(self->ndim(), tpl));
                 return py::make_tuple(std::get<0>(splitted),
                                       shape_tuple(std::get<1>(splitted)),
                                       shape_tuple(std::get<2>(splitted)));
               },
+              py::arg("tpl"),
               "Split an index into a spin block descriptor, a spatial block index and an "
               "in-block index.")
         .def(
               "combine",
-              [](std::shared_ptr<MoIndexTranslation> self, py::tuple bidx,
-                 py::tuple ibidx) {
+              [](std::shared_ptr<MoIndexTranslation> self, IdxTuple bidx,
+                 IdxTuple ibidx) {
                 return shape_tuple(self->combine(parse_tuple(self->ndim(), bidx),
                                                  parse_tuple(self->ndim(), ibidx)));
               },
+              py::arg("bidx"), py::arg("ibidx"),
               "Combine a block index and an in-block index into the appropriate index. "
               "Effectively undoes the effect of 'split'.")
         .def(
               "combine",
               [](std::shared_ptr<MoIndexTranslation> self, std::string spin_block,
-                 py::tuple bidx, py::tuple ibidx) {
+                 IdxTuple bidx, IdxTuple ibidx) {
                 return shape_tuple(self->combine(spin_block,
                                                  parse_tuple(self->ndim(), bidx),
                                                  parse_tuple(self->ndim(), ibidx)));
               },
+              py::arg("spin_block"), py::arg("bidx"), py::arg("ibidx"),
               "Combine a spin block (given as a string of 'a's or 'b's), a spatial-only "
               "block index and an in-block index into the appropriate index. Essentially "
               "undoes the effect of 'spin_of', 'block_index_spatial_of' and "
               "'inblock_index_of'.")
         .def(
               "hf_provider_index_of",
-              [](std::shared_ptr<MoIndexTranslation> self, py::tuple index) {
+              [](std::shared_ptr<MoIndexTranslation> self, IdxTuple index) {
                 return shape_tuple(
                       self->hf_provider_index_of(parse_tuple(self->ndim(), index)));
               },
+              py::arg("index"),
               "Map an index (given in the space passed upon construction) to the "
               "indexing "
               "convention of the host program provided to adcc as the HF provider.")
         //
         .def("map_range_to_hf_provider", &MoIndexTranslation_map_range_to_hf_provider,
+             py::arg("ranges"),
              "Map a range of indices to host program indices, i.e. the indexing "
              "convention\n"
              "used in the HfProvider, which provides the SCF data to adcc.\n"
              "\n"
              "Since the mapping between subspace and host program indices might not be "
              "contiguous,\n"
-             "a list of pairs of ranges is returned. In each pair, the first entry "
-             "represents a\n"
-             "range of indices (indexed in the MO subspace) and the second entry "
-             "represents\n"
-             "the equivalent range of indices in the Hartree-Fock provider these are "
-             "mapped to.\n"
+             "a list of range mappings is returned. Each mapping is a dict with the "
+             "keys\n"
+             "'from' and 'to': 'from' holds a range of indices (indexed in the MO "
+             "subspace)\n"
+             "and 'to' the equivalent range of indices in the Hartree-Fock provider "
+             "these are\n"
+             "mapped to. Both are given as a tuple of half-open [start, end) index "
+             "pairs, one\n"
+             "pair for each dimension, i.e. in the same format as the `ranges` "
+             "argument.\n"
+             "\n"
+             "For example, for a two-dimensional space:\n"
+             "  [{'from': ((0, 1), (0, 1)), 'to': ((0, 1), (2, 3))}, ...]\n"
              "\n"
              "  ranges    Tuple of pairs of indices: One index pair for each dimension. "
              "Each\n"
diff --git a/libadcc_src/pyiface/export_MoSpaces.cc b/libadcc_src/pyiface/export_MoSpaces.cc
index a301a2ee..a8781341 100644
--- a/libadcc_src/pyiface/export_MoSpaces.cc
+++ b/libadcc_src/pyiface/export_MoSpaces.cc
@@ -32,10 +32,12 @@ void export_MoSpaces(py::module& m) {
         "information about them. Python binding to :cpp:class:`libadcc::MoSpaces`.")
         .def(py::init<const HartreeFockSolution_i&, std::shared_ptr<const AdcMemory>,
                       std::vector<size_t>, std::vector<size_t>, std::vector<size_t>>(),
+             py::arg("hf"), py::arg("adcmem"), py::arg("core_orbitals"),
+             py::arg("frozen_core_orbitals"), py::arg("frozen_virtuals"),
              "Construct an MoSpaces object from a HartreeFockSolution_i, a pointer to\n"
              "an AdcMemory object.\n"
              "\n"
-             "adcmem_ptr        ADC memory keep-alive object to be used in all Tensors\n"
+             "adcmem            ADC memory keep-alive object to be used in all Tensors\n"
              "                  constructed using this MoSpaces object.\n"
              "core_orbitals     List of orbitals indices (in the full fock space, "
              "original\n"
@@ -56,11 +58,11 @@ void export_MoSpaces(py::module& m) {
              "beta\n"
              "                  orbitals has to be selected.\n")
         //
-        .def("n_orbs", &MoSpaces::n_orbs,
+        .def("n_orbs", &MoSpaces::n_orbs, py::arg("space"),
              "The number of orbitals in a particular orbital subspace")
-        .def("n_orbs_alpha", &MoSpaces::n_orbs_alpha,
+        .def("n_orbs_alpha", &MoSpaces::n_orbs_alpha, py::arg("space"),
              "The number of alpha orbitals in a particular orbital subspace")
-        .def("n_orbs_beta", &MoSpaces::n_orbs_beta,
+        .def("n_orbs_beta", &MoSpaces::n_orbs_beta, py::arg("space"),
              "The number of beta orbitals in a particular orbital subspace")
         //
         .def_readonly("point_group", &MoSpaces::point_group,
diff --git a/libadcc_src/pyiface/export_ReferenceState.cc b/libadcc_src/pyiface/export_ReferenceState.cc
index 2a34fb2c..907875c7 100644
--- a/libadcc_src/pyiface/export_ReferenceState.cc
+++ b/libadcc_src/pyiface/export_ReferenceState.cc
@@ -19,6 +19,7 @@
 
 #include "../ReferenceState.hh"
 #include "hartree_fock_solution_hack.hh"
+#include "ndarray.hh"
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
@@ -55,19 +56,21 @@ py::object convert_timer(const Timer& timer) {
 
 void export_ReferenceState(py::module& m) {
 
-  py::class_<ReferenceState, std::shared_ptr<ReferenceState>>(
+  py::class_<ReferenceState, std::shared_ptr<ReferenceState>> reference_state(
         m, "ReferenceState",
         "Class representing information about the reference state for adcc. Python "
         "binding to"
-        ":cpp:class:`libadcc::ReferenceState`.")
+        ":cpp:class:`libadcc::ReferenceState`.");
+  reference_state
         .def(py::init<std::shared_ptr<const HartreeFockSolution_i>,
                       std::shared_ptr<const MoSpaces>, bool>(),
+             py::arg("hfsoln"), py::arg("mo"), py::arg("symmetry_check_on_import"),
              "Setup a ReferenceStateject using an MoSpaces object.\n"
              "\n"
-             "hfsoln_ptr        Pointer to the Interface to the host program,\n"
+             "hfsoln            Pointer to the Interface to the host program,\n"
              "                  providing the HartreeFockSolution data, which\n"
              "                  will be provided by this object.\n"
-             "mo_ptr            MoSpaces object containing info about the MoSpace setup\n"
+             "mo                MoSpaces object containing info about the MoSpace setup\n"
              "                  and the point group symmetry.\n"
              "symmetry_check_on_import\n"
              "                  Should symmetry of the imported objects be checked\n"
@@ -112,25 +115,30 @@ void export_ReferenceState(py::module& m) {
               "nuclear_total_charge",
               [](const ReferenceState& ref) { return ref.nuclear_multipole(0)[0]; })
         .def_property_readonly("nuclear_dipole",
-                               [](const ReferenceState& ref) {
+                               [](const ReferenceState& ref) -> NDArray<scalar_type, 1> {
                                  py::array_t<scalar_type> ret(std::vector<ssize_t>{3});
                                  auto res = ref.nuclear_multipole(1);
                                  std::copy(res.begin(), res.end(), ret.mutable_data());
                                  return ret;
                                })
-        .def("nuclear_quadrupole",
-             [](const ReferenceState& ref, std::array<scalar_type, 3> gauge_origin) {
-               py::array_t<scalar_type> ret(std::vector<ssize_t>{6});
-               auto res = ref.nuclear_multipole(2, gauge_origin);
-               std::copy(res.begin(), res.end(), ret.mutable_data());
-               return ret;
-             })
-        .def("gauge_origin_to_xyz",
-             [](const ReferenceState& ref, std::string gauge_origin) {
-               auto vec = ref.gauge_origin_to_xyz(gauge_origin);
-               // Make sure a tuple is returned
-               return py::make_tuple(vec[0], vec[1], vec[2]);
-             })
+        .def(
+              "nuclear_quadrupole",
+              [](const ReferenceState& ref,
+                 std::tuple<scalar_type, scalar_type, scalar_type> gauge_origin)
+                    -> NDArray<scalar_type, 1> {
+                py::array_t<scalar_type> ret(std::vector<ssize_t>{6});
+                auto res = ref.nuclear_multipole(2, gauge_origin);
+                std::copy(res.begin(), res.end(), ret.mutable_data());
+                return ret;
+              },
+              py::arg("gauge_origin"))
+        .def(
+              "gauge_origin_to_xyz",
+              [](const ReferenceState& ref, std::string gauge_origin)
+                    -> std::tuple<scalar_type, scalar_type, scalar_type> {
+                return ref.gauge_origin_to_xyz(gauge_origin);
+              },
+              py::arg("gauge_origin"))
         .def_property_readonly("conv_tol", &ReferenceState::conv_tol,
                                "SCF convergence tolererance")
         .def_property_readonly("energy_scf", &ReferenceState::energy_scf,
@@ -139,20 +147,23 @@ void export_ReferenceState(py::module& m) {
                                &ReferenceState::nuclear_repulsion_energy,
                                "The nuclear repulsion energy")
         //
-        .def("orbital_energies", &ReferenceState::orbital_energies,
+        .def("orbital_energies", &ReferenceState::orbital_energies, py::arg("space"),
              "Return the orbital energies corresponding to the provided space")
         .def("orbital_coefficients", &ReferenceState::orbital_coefficients,
+             py::arg("space"),
              "Return the molecular orbital coefficients corresponding to the provided "
              "space (alpha and beta coefficients are returned)")
         .def("orbital_coefficients_alpha", &ReferenceState::orbital_coefficients_alpha,
+             py::arg("space"),
              "Return the alpha molecular orbital coefficients corresponding to the "
              "provided space")
         .def("orbital_coefficients_beta", &ReferenceState::orbital_coefficients_beta,
+             py::arg("space"),
              "Return the beta molecular orbital coefficients corresponding to the "
              "provided space")
-        .def("fock", &ReferenceState::fock,
+        .def("fock", &ReferenceState::fock, py::arg("space"),
              "Return the Fock matrix block corresponding to the provided space.")
-        .def("eri", &ReferenceState::eri,
+        .def("eri", &ReferenceState::eri, py::arg("space"),
              "Return the ERI (electron-repulsion integrals) tensor block corresponding "
              "to the provided space.")
         //
@@ -162,18 +173,22 @@ void export_ReferenceState(py::module& m) {
              "is requested by a call to above fock() or eri() functions. This function "
              "call, however, instructs the class to immediately import *all* such "
              "blocks. Typically you do not want to do this.")
-        .def_property("cached_fock_blocks", &ReferenceState::cached_fock_blocks,
-                      &ReferenceState::set_cached_fock_blocks,
-                      "Get or set the list of momentarily cached Fock matrix blocks\n"
-                      "\n"
-                      "Setting this property allows to drop fock matrix blocks if they "
-                      "are no longer needed to save memory.")
-        .def_property("cached_eri_blocks", &ReferenceState::cached_eri_blocks,
-                      &ReferenceState::set_cached_eri_blocks,
-                      "Get or set the list of momentarily cached ERI tensor blocks\n"
-                      "\n"
-                      "Setting this property allows to drop ERI tensor blocks if they "
-                      "are no longer needed to save memory.")
+        .def_property(
+              "cached_fock_blocks", &ReferenceState::cached_fock_blocks,
+              py::cpp_function(&ReferenceState::set_cached_fock_blocks,
+                               py::is_method(reference_state), py::arg("newlist")),
+              "Get or set the list of momentarily cached Fock matrix blocks\n"
+              "\n"
+              "Setting this property allows to drop fock matrix blocks if they "
+              "are no longer needed to save memory.")
+        .def_property(
+              "cached_eri_blocks", &ReferenceState::cached_eri_blocks,
+              py::cpp_function(&ReferenceState::set_cached_eri_blocks,
+                               py::is_method(reference_state), py::arg("newlist")),
+              "Get or set the list of momentarily cached ERI tensor blocks\n"
+              "\n"
+              "Setting this property allows to drop ERI tensor blocks if they "
+              "are no longer needed to save memory.")
         .def("flush_hf_cache", &ReferenceState::flush_hf_cache,
              "Tell the contained HartreeFockSolution_i object (which was passed upon "
              "construction), that a larger amount of import operations is done and that "
diff --git a/libadcc_src/pyiface/export_Symmetry.cc b/libadcc_src/pyiface/export_Symmetry.cc
index 05e3dff5..02d9abd3 100644
--- a/libadcc_src/pyiface/export_Symmetry.cc
+++ b/libadcc_src/pyiface/export_Symmetry.cc
@@ -29,14 +29,17 @@ namespace py = pybind11;
 
 void export_Symmetry(py::module& m) {
 
-  py::class_<Symmetry, std::shared_ptr<Symmetry>>(
-        m, "Symmetry", "Container for Tensor symmetry information")
+  py::class_<Symmetry, std::shared_ptr<Symmetry>> symmetry(
+        m, "Symmetry", "Container for Tensor symmetry information");
+  symmetry
         .def(py::init<std::shared_ptr<const MoSpaces>, const std::string&>(),
+             py::arg("mospaces"), py::arg("space"),
              "Construct a Symmetry class from an MoSpaces object and the identifier for "
              "the space (e.g. o1o1, v1o1, o3v2o1v1, ...). Python binding to "
              ":cpp:class:`libadcc::Symmetry`.")
         .def(py::init<std::shared_ptr<const MoSpaces>, const std::string&,
                       std::map<std::string, std::pair<size_t, size_t>>>(),
+             py::arg("mospaces"), py::arg("space"), py::arg("extra_axes_orbs"),
              "Construct a Symmetry class from an MoSpaces object, a space string and a "
              "map to supply the number of orbitals for some additional axes.\nFor the "
              "additional axis the pair contains either two numbers (for the number of "
@@ -60,12 +63,15 @@ void export_Symmetry(py::module& m) {
         .def("describe", &Symmetry::describe, "Return a descriptive string.")
         //
         .def_property("irreps_allowed", &Symmetry::irreps_allowed,
-                      &Symmetry::set_irreps_allowed,
+                      py::cpp_function(&Symmetry::set_irreps_allowed,
+                                       py::is_method(symmetry), py::arg("irreps")),
                       "The list of irreducible representations, for which the tensor "
                       "shall be non-zero. If this is *not* set, i.e. an empty list, all "
                       "irreps will be allowed.")
         .def_property(
-              "permutations", &Symmetry::permutations, &Symmetry::set_permutations,
+              "permutations", &Symmetry::permutations,
+              py::cpp_function(&Symmetry::set_permutations, py::is_method(symmetry),
+                               py::arg("permutations")),
               "The list of index permutations, which do not change the tensor.\n"
               "A minus may be used to indicate anti-symmetric\n"
               "permutations with respect to the first (reference) permutation.\n"
@@ -76,13 +82,15 @@ void export_Symmetry(py::module& m) {
               "the symmetry. Beware that the check for errors and conflicts\n"
               "is only rudimentary at the moment.")
         .def_property("spin_block_maps", &Symmetry::spin_block_maps,
-                      &Symmetry::set_spin_block_maps,
+                      py::cpp_function(&Symmetry::set_spin_block_maps,
+                                       py::is_method(symmetry), py::arg("spin_maps")),
                       "A list of tuples of the form (\"aaaa\", \"bbbb\", -1.0), i.e.\n"
                       "two spin blocks followed by a factor. This maps the second onto "
                       "the first\n"
                       "with a factor of -1.0 between them.")
         .def_property("spin_blocks_forbidden", &Symmetry::spin_blocks_forbidden,
-                      &Symmetry::set_spin_blocks_forbidden,
+                      py::cpp_function(&Symmetry::set_spin_blocks_forbidden,
+                                       py::is_method(symmetry), py::arg("forbidden")),
                       "List of spin-blocks, which are marked forbidden (i.e. enforce "
                       "them to stay zero).\n"
                       "Blocks are given as a string in the letters 'a' and 'b', e.g. "
@@ -94,6 +102,7 @@ void export_Symmetry(py::module& m) {
   // Factories for common cases
   //
   m.def("make_symmetry_orbital_energies", &make_symmetry_orbital_energies,
+        py::arg("mospaces"), py::arg("space"),
         "Return the Symmetry object like it would be set up for the passed subspace \n"
         "of the orbital energies tensor.\n"
         "\n"
@@ -101,6 +110,7 @@ void export_Symmetry(py::module& m) {
         "  space       space string (e.g. o1)");
 
   m.def("make_symmetry_orbital_coefficients", &make_symmetry_orbital_coefficients,
+        py::arg("mospaces"), py::arg("space"), py::arg("n_bas"), py::arg("blocks") = "ab",
         "Return the Symmetry object like it would be set up for the passed subspace \n"
         "of the orbital coefficients tensor.\n"
         "\n"
@@ -109,25 +119,28 @@ void export_Symmetry(py::module& m) {
         "  n_bas       Number of basis functions\n"
         "  blocks      Spin blocks to include. Valid are \"ab\", \"a\" and \"b\".");
 
-  m.def("make_symmetry_eri", &make_symmetry_eri,
+  m.def("make_symmetry_eri", &make_symmetry_eri, py::arg("mospaces"), py::arg("space"),
         "Return the Symmetry object like it would be set up for the passed subspace \n"
         "of the electron-repulsion tensor.\n"
         "\n"
         "  mospaces    MoSpaces object\n"
         "  space       Space string (e.g. o1v1o1v1)\n");
 
-  m.def("make_symmetry_operator", &make_symmetry_operator,
+  m.def("make_symmetry_operator", &make_symmetry_operator, py::arg("mospaces"),
+        py::arg("space"), py::arg("operator_symmetry"),
+        py::arg("cartesian_transformation"),
         "Return the Symmetry object for an orbital subspace block of a one-particle "
         "operator\n"
         "\n"
         "  mospaces    MoSpaces object\n"
         "  space       Space string (e.g. o1v1)\n"
-        "  symmetry    Describes the symmetry of the tensor (only in effect if both \n"
+        "  operator_symmetry\n"
+        "              Describes the symmetry of the tensor (only in effect if both \n"
         "              subspaces of the space string are identical).\n"
         "              Valid are \"nosymmetry\", \"hermitian\" and \"antihermitian\".\n"
         "  cartesian_transformation\n"
         "              The cartesian function according to which the operator "
-        "transforms.\n"
+        "              transforms.\n"
         "\n"
         "Valid cartesian_transformation values include:\n"
         "     \"1\"                   Totally symmetric (default)\n"
@@ -136,6 +149,8 @@ void export_Symmetry(py::module& m) {
         "     \"Rx\", \"Ry\", \"Rz\"      Rotations about the coordinate axis\n");
 
   m.def("make_symmetry_operator_basis", &make_symmetry_operator_basis,
+        py::arg("mospaces"), py::arg("n_bas"), py::arg("operator_symmetry"),
+        py::arg("n_particle_op"), py::arg("blocks"),
         "Return the symmetry object for an operator in the AO basis. The object will\n"
         "represent a block-diagonal matrix of the form\n"
         "    ( M 0 )\n"
@@ -143,17 +158,18 @@ void export_Symmetry(py::module& m) {
         "where M is an n_bas x n_bas block and is indentical in upper-left\n"
         "and lower-right.\n"
         "\n"
-        "mospaces_ptr      MoSpaces pointer\n"
+        "mospaces          MoSpaces pointer\n"
         "n_bas             Number of AO basis functions\n"
         "operator_symmetry Is the tensor symmetric (hermitian/antihermitian, only\n"
         "                  in effect if both space axes identical).\n"
-        "                  Nosymmetry disables a setup of permutational symmetry.\n"
+        "                  Nosymmetry disables a setup of 'bra-ket' symmetry.\n"
         "n_particle_op     NParticleOperator\n"
         "blocks            Which blocks of the operator to return. Valid values\n"
         "                  are 'ab' to return a tensor for both alpha and beta\n"
         "                  block as a block-diagonal tensor, 'a' to only return a\n"
         "                  tensor for only alpha block.\n");
-  m.def("make_symmetry_triples", &make_symmetry_triples,
+  m.def("make_symmetry_triples", &make_symmetry_triples, py::arg("mospaces"),
+        py::arg("space"),
         "Return the Symmetry object like it would be set up for the passed subspace \n"
         "of a triples amplitude tensor.\n"
         "\n"
diff --git a/libadcc_src/pyiface/export_Tensor.cc b/libadcc_src/pyiface/export_Tensor.cc
index ecff8ac7..bb32172e 100644
--- a/libadcc_src/pyiface/export_Tensor.cc
+++ b/libadcc_src/pyiface/export_Tensor.cc
@@ -19,19 +19,24 @@
 
 #include "../Tensor.hh"
 #include "../exceptions.hh"
+#include "ndarray.hh"
 #include "util.hh"
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
+#include <pybind11/typing.h>
 #include <sstream>
 
 namespace libadcc {
 
 namespace py = pybind11;
-using namespace pybind11::literals;
 typedef std::shared_ptr<Tensor> ten_ptr;
 
+using Permutations =
+      py::typing::Union<py::typing::Iterable<py::int_>,
+                        py::typing::Iterable<py::typing::Iterable<py::int_>>>;
+
 static std::vector<std::vector<size_t>> parse_permutations(
-      const py::iterable& permutations) {
+      const Permutations& permutations) {
   bool iterator_of_ints = true;
   for (auto tpl : permutations) {
     if (!py::isinstance<py::int_>(tpl)) {
@@ -103,7 +108,7 @@ static std::vector<size_t> convert_index_tuple(const ten_ptr& self, py::tuple id
 //
 // Extra defs
 //
-static py::tuple Tensor_shape(const Tensor& self) { return shape_tuple(self.shape()); }
+static auto Tensor_shape(const Tensor& self) { return shape_tuple(self.shape()); }
 
 static py::array_t<scalar_type> Tensor_to_ndarray(const Tensor& self) {
   // Get an empty array of the required shape and export the data into it.
@@ -113,8 +118,7 @@ static py::array_t<scalar_type> Tensor_to_ndarray(const Tensor& self) {
 }
 
 static ten_ptr Tensor_from_ndarray_tol(
-      ten_ptr self,
-      py::array_t<scalar_type, py::array::c_style | py::array::forcecast> in_array,
+      ten_ptr self, py::array_t<scalar_type, py::array::c_style> in_array,
       double symmetry_tolerance) {
   py::ssize_t nd = in_array.ndim();
   if (nd < 1) throw invalid_argument("Cannot import from 0D array.");
@@ -127,7 +131,8 @@ static ten_ptr Tensor_from_ndarray_tol(
   return self;
 }
 
-static ten_ptr Tensor_from_ndarray(ten_ptr self, py::array in_array) {
+static ten_ptr Tensor_from_ndarray(
+      ten_ptr self, py::array_t<scalar_type, py::array::c_style> in_array) {
   return Tensor_from_ndarray_tol(self, in_array, 0.0);
 }
 
@@ -140,8 +145,9 @@ static scalar_type Tensor_dot(const Tensor& self, ten_ptr other) {
   return self.dot({other})[0];
 }
 
-static py::array_t<scalar_type> Tensor_dot_list(const Tensor& self, py::list tensors) {
-  std::vector<ten_ptr> parsed   = extract_tensors(tensors);
+static NDArray<scalar_type, 1> Tensor_dot_list(const Tensor& self,
+                                               py::typing::List<Tensor> tensors) {
+  std::vector<ten_ptr> parsed   = extract_tensors<py::list>(tensors);
   std::vector<scalar_type> dots = self.dot(parsed);
   py::array_t<scalar_type> ret(dots.size());
   std::copy(dots.begin(), dots.end(), ret.mutable_data());
@@ -157,7 +163,8 @@ static ten_ptr Tensor_transpose_1(const Tensor& self) {
   return self.transpose(vec_axes);
 }
 
-static ten_ptr Tensor_transpose_2(const Tensor& self, py::tuple axes) {
+static ten_ptr Tensor_transpose_2(const Tensor& self,
+                                  py::typing::Tuple<py::ssize_t, py::ellipsis> axes) {
   std::vector<size_t> vec_axes(py::len(axes));
   for (size_t i = 0; i < py::len(axes); ++i) {
     vec_axes[i] = axes[i].cast<size_t>();
@@ -165,11 +172,11 @@ static ten_ptr Tensor_transpose_2(const Tensor& self, py::tuple axes) {
   return self.transpose(vec_axes);
 }
 
-static ten_ptr Tensor_symmetrise_1(const Tensor& self, py::list permutations) {
+static ten_ptr Tensor_symmetrise_1(const Tensor& self, Permutations permutations) {
   return self.symmetrise(parse_permutations(permutations));
 }
 
-static ten_ptr Tensor_symmetrise_2(const Tensor& self, py::args permutations) {
+static ten_ptr Tensor_symmetrise_2(const Tensor& self, py::Args<py::int_> permutations) {
   if (py::len(permutations) == 0) {
     if (self.ndim() != 2) {
       throw invalid_argument(
@@ -180,11 +187,12 @@ static ten_ptr Tensor_symmetrise_2(const Tensor& self, py::args permutations) {
   return self.symmetrise(parse_permutations(permutations));
 }
 
-static ten_ptr Tensor_antisymmetrise_1(const Tensor& self, py::list permutations) {
+static ten_ptr Tensor_antisymmetrise_1(const Tensor& self, Permutations permutations) {
   return self.antisymmetrise(parse_permutations(permutations));
 }
 
-static ten_ptr Tensor_antisymmetrise_2(const Tensor& self, py::args permutations) {
+static ten_ptr Tensor_antisymmetrise_2(const Tensor& self,
+                                       py::Args<py::int_> permutations) {
   if (py::len(permutations) == 0) {
     if (self.ndim() != 2) {
       throw invalid_argument(
@@ -195,7 +203,8 @@ static ten_ptr Tensor_antisymmetrise_2(const Tensor& self, py::args permutations
   return self.antisymmetrise(parse_permutations(permutations));
 }
 
-static py::object tensordot_1(ten_ptr a, ten_ptr b, py::iterable axes) {
+static py::typing::Union<ten_ptr, scalar_type> tensordot_1(
+      ten_ptr a, ten_ptr b, py::typing::Iterable<py::typing::Iterable<py::int_>> axes) {
   if (py::len(axes) != 2) {
     throw invalid_argument("axes needs to be an iterable of length 2");
   }
@@ -215,7 +224,8 @@ static py::object tensordot_1(ten_ptr a, ten_ptr b, py::iterable axes) {
     return py::cast(res.tensor_ptr);
   }
 }
-static py::object tensordot_2(ten_ptr a, ten_ptr b, size_t axes) {
+static py::typing::Union<ten_ptr, scalar_type> tensordot_2(ten_ptr a, ten_ptr b,
+                                                           size_t axes) {
   std::vector<size_t> a_axes;
   std::vector<size_t> b_axes;
   for (size_t i = 0; i < axes; ++i) {
@@ -231,9 +241,11 @@ static py::object tensordot_2(ten_ptr a, ten_ptr b, size_t axes) {
   }
 }
 
-static py::object tensordot_3(ten_ptr a, ten_ptr b) { return tensordot_2(a, b, 2); }
+static py::typing::Union<ten_ptr, scalar_type> tensordot_3(ten_ptr a, ten_ptr b) {
+  return tensordot_2(a, b, 2);
+}
 
-static ten_ptr Tensor_diagonal(ten_ptr ten, py::args permutations) {
+static ten_ptr Tensor_diagonal(ten_ptr ten, py::Args<py::int_> permutations) {
   std::vector<size_t> axes;
   if (py::len(permutations) == 0) {
     axes.push_back(0);
@@ -244,12 +256,12 @@ static ten_ptr Tensor_diagonal(ten_ptr ten, py::args permutations) {
   return ten->diagonal(axes);
 }
 
-static ten_ptr direct_sum(ten_ptr a, ten_ptr b) { return a->direct_sum(b); }
+static auto direct_sum(ten_ptr a, ten_ptr b) { return a->direct_sum(b); }
 
-static double Tensor_trace_1(std::string subscripts, const Tensor& tensor) {
+static auto Tensor_trace_1(std::string subscripts, const Tensor& tensor) {
   return tensor.trace(subscripts);
 }
-static double Tensor_trace_2(const Tensor& tensor) {
+static auto Tensor_trace_2(const Tensor& tensor) {
   if (tensor.ndim() != 2) {
     throw invalid_argument(
           "trace function without arguments may only be used for matrices.");
@@ -258,7 +270,8 @@ static double Tensor_trace_2(const Tensor& tensor) {
 }
 
 static ten_ptr linear_combination_strict(
-      py::array_t<scalar_type, py::array::c_style> coefficients, py::list tensors) {
+      py::array_t<scalar_type, py::array::c_style> coefficients,
+      py::typing::List<Tensor> tensors) {
 
   if (coefficients.ndim() != 1) {
     throw invalid_argument("coefficients array needs to have exactly one dimension.");
@@ -267,7 +280,9 @@ static ten_ptr linear_combination_strict(
   const scalar_type* in_data = coefficients.data();
   std::vector<scalar_type> scalars(in_size);
   std::copy(in_data, in_data + in_size, scalars.data());
-  std::vector<ten_ptr> parsed = extract_tensors(tensors);
+  std::vector<ten_ptr> parsed = extract_tensors<py::list>(tensors);
+
+  if (parsed.empty()) throw runtime_error("No tensors parsed during linear combination");
 
   auto ret = parsed[0]->zeros_like();
   ret->add_linear_combination(scalars, parsed);
@@ -278,46 +293,51 @@ static ten_ptr linear_combination_strict(
 // Element access
 //
 
-static py::list Tensor_select_n_min(const ten_ptr& self, size_t n) {
+using ElementList =
+      py::typing::List<py::typing::Tuple<py::typing::List<size_t>, scalar_type>>;
+
+static ElementList Tensor_select_n_min(const ten_ptr& self, size_t n) {
   std::vector<std::pair<std::vector<size_t>, scalar_type>> ret = self->select_n_min(n);
   py::list li;
   for (auto p : ret) li.append(py::make_tuple(p.first, p.second));
   return li;
 }
 
-static py::list Tensor_select_n_max(const ten_ptr& self, size_t n) {
+static ElementList Tensor_select_n_max(const ten_ptr& self, size_t n) {
   std::vector<std::pair<std::vector<size_t>, scalar_type>> ret = self->select_n_max(n);
   py::list li;
   for (auto p : ret) li.append(py::make_tuple(p.first, p.second));
   return li;
 }
 
-static py::list Tensor_select_n_absmin(const ten_ptr& self, size_t n) {
+static ElementList Tensor_select_n_absmin(const ten_ptr& self, size_t n) {
   std::vector<std::pair<std::vector<size_t>, scalar_type>> ret = self->select_n_absmin(n);
   py::list li;
   for (auto p : ret) li.append(py::make_tuple(p.first, p.second));
   return li;
 }
 
-static py::list Tensor_select_n_absmax(const ten_ptr& self, size_t n) {
+static ElementList Tensor_select_n_absmax(const ten_ptr& self, size_t n) {
   std::vector<std::pair<std::vector<size_t>, scalar_type>> ret = self->select_n_absmax(n);
   py::list li;
   for (auto p : ret) li.append(py::make_tuple(p.first, p.second));
   return li;
 }
 
-static bool Tensor_is_allowed(const ten_ptr& self, py::tuple idcs) {
+static bool Tensor_is_allowed(const ten_ptr& self,
+                              py::typing::Tuple<py::ssize_t, py::ellipsis> idcs) {
   return self->is_element_allowed(convert_index_tuple(self, idcs));
 }
 
-static scalar_type Tensor__getitem__(const ten_ptr& self, py::tuple idcs) {
+static scalar_type Tensor__getitem__(const ten_ptr& self,
+                                     py::typing::Tuple<py::ssize_t, py::ellipsis> idcs) {
   return self->get_element(convert_index_tuple(self, idcs));
 }
 
-static scalar_type Tensor__setitem__(const ten_ptr& self, py::tuple idcs,
-                                     scalar_type value) {
+static void Tensor__setitem__(const ten_ptr& self,
+                              py::typing::Tuple<py::ssize_t, py::ellipsis> idcs,
+                              scalar_type value) {
   self->set_element(convert_index_tuple(self, idcs), value);
-  return value;
 }
 
 //
@@ -325,7 +345,7 @@ static scalar_type Tensor__setitem__(const ten_ptr& self, py::tuple idcs,
 //    See https://docs.python.org/3/library/operator.html for details
 //
 
-static py::object Tensor___str__(const Tensor& self) {
+static py::str Tensor___str__(const Tensor& self) {
   if (self.size() < 50000) {
     return py::str(Tensor_to_ndarray(self));
   } else {
@@ -340,7 +360,7 @@ static py::object Tensor___str__(const Tensor& self) {
   }
 }
 
-static py::object Tensor___repr__(const Tensor& self) {
+static py::str Tensor___repr__(const Tensor& self) {
   // TODO extremely rudimentary information for now
   //      goal would be an unambiguous representation instead
   return Tensor___str__(self);
@@ -422,10 +442,10 @@ static ten_ptr Tensor__matmul__(const ten_ptr& self, const ten_ptr& other) {
 }
 
 void export_Tensor(py::module& m) {
-  py::class_<Tensor, std::shared_ptr<Tensor>>(
+  py::class_<Tensor, std::shared_ptr<Tensor>> tensor(
         m, "Tensor",
-        "Class representing the Tensor objects used for computations in adcc")
-        .def(py::init(&make_tensor_zero),
+        "Class representing the Tensor objects used for computations in adcc");
+  tensor.def(py::init(&make_tensor_zero), py::arg("symmetry"),
              "Construct a Tensor object using a Symmetry object describing its symmetry "
              "properties.\n"
              "The returned object is not guaranteed to contain initialised memory. "
@@ -435,7 +455,9 @@ void export_Tensor(py::module& m) {
         .def_property_readonly("size", &Tensor::size)
         .def_property_readonly("space", &Tensor::space)
         .def_property_readonly("subspaces", &Tensor::subspaces)
-        .def_property("flags", &Tensor::flags, &Tensor::set_flags)
+        .def_property("flags", &Tensor::flags,
+                      py::cpp_function(&Tensor::set_flags, py::is_method(tensor),
+                                       py::arg("new_flags")))
         //
         .def_property_readonly("needs_evaluation", &Tensor::needs_evaluation,
                                "Does the tensor need evaluation or is it fully evaluated "
@@ -455,52 +477,53 @@ void export_Tensor(py::module& m) {
         .def("set_random", &Tensor_set_random,
              "Set all tensor elements to random data, adhering to the internal "
              "symmetry.")
-        .def("set_mask", &Tensor::set_mask,
+        .def("set_mask", &Tensor::set_mask, py::arg("mask"), py::arg("value"),
              "Set all elements corresponding to an index mask, which is given by a "
              "string eg. 'iijkli' sets elements T_{iijkli}")
         .def("diagonal", &Tensor_diagonal)
         .def("copy", &Tensor::copy, "Returns a deep copy of the tensor.")
-        .def("dot", &Tensor_dot)
-        .def("dot", &Tensor_dot_list)
+        .def("dot", &Tensor_dot, py::arg("other"))
+        .def("dot", &Tensor_dot_list, py::arg("tensors"))
         .def_property_readonly("T", &Tensor_transpose_1)
         .def("transpose", &Tensor_transpose_1)
-        .def("transpose", &Tensor_transpose_2)
-        .def("symmetrise", &Tensor_symmetrise_1)
+        .def("transpose", &Tensor_transpose_2, py::arg("axes"))
+        .def("symmetrise", &Tensor_symmetrise_1, py::arg("permutations"))
         .def("symmetrise", &Tensor_symmetrise_2)
-        .def("antisymmetrise", &Tensor_antisymmetrise_1)
+        .def("antisymmetrise", &Tensor_antisymmetrise_1, py::arg("permutations"))
         .def("antisymmetrise", &Tensor_antisymmetrise_2)
         .def("to_ndarray", &Tensor_to_ndarray,
              "Export the tensor data to a standard np::ndarray by making a copy.")
-        .def("set_from_ndarray", &Tensor_from_ndarray,
-             "Set all tensor elements from a standard np::ndarray by making a copy. "
-             "Provide an optional tolerance argument to increase the tolerance for the "
-             "check for symmetry consistency.")
-        .def("set_from_ndarray", &Tensor_from_ndarray_tol,
+        .def("set_from_ndarray", &Tensor_from_ndarray, py::arg("in_array"),
+             "Set all tensor elements from a standard np::ndarray by making a copy.")
+        .def("set_from_ndarray", &Tensor_from_ndarray_tol, py::arg("in_array"),
+             py::arg("symmetry_tolerance"),
              "Set all tensor elements from a standard np::ndarray by making a copy. "
              "Provide an optional tolerance argument to increase the tolerance for the "
              "check for symmetry consistency.")
         .def("describe_symmetry", &Tensor::describe_symmetry,
              "Return a string providing a hopefully descriptive representation of the "
              "symmetry information stored inside the tensor.")
-        .def("describe_expression", &Tensor::describe_expression,
+        .def("describe_expression", &Tensor::describe_expression, py::arg("stage"),
              "Return a string providing a hopefully descriptive representation of the "
              "tensor expression stored inside the object.")
         .def("describe_expression",
              [](ten_ptr t) { return t->describe_expression("unoptimised"); })
         //
-        .def("__getitem__", &Tensor__getitem__,
+        .def("__getitem__", &Tensor__getitem__, py::arg("idcs"),
              "Get a tensor element or a slice of tensor elements.")
-        .def("__setitem__", &Tensor__setitem__,
+        .def("__setitem__", &Tensor__setitem__, py::arg("idcs"), py::arg("value"),
              "Set a tensor element or a slice of tensor elements. The operation will "
              "adhere symmetry, i.e. alter all elements equivalent by symmetry at once.")
-        .def("is_allowed", &Tensor_is_allowed,
+        .def("is_allowed", &Tensor_is_allowed, py::arg("idcs"),
              " Is a particular index allowed by symmetry")
-        .def("select_n_absmax", &Tensor_select_n_absmax,
+        .def("select_n_absmax", &Tensor_select_n_absmax, py::arg("n"),
              "Select the n absolute maximal elements.")
-        .def("select_n_absmin", &Tensor_select_n_absmin,
+        .def("select_n_absmin", &Tensor_select_n_absmin, py::arg("n"),
              "Select the n absolute minimal elements.")
-        .def("select_n_max", &Tensor_select_n_max, "Select the n maximal elements.")
-        .def("select_n_min", &Tensor_select_n_min, "Select the n minimal elements.")
+        .def("select_n_max", &Tensor_select_n_max, py::arg("n"),
+             "Select the n maximal elements.")
+        .def("select_n_min", &Tensor_select_n_min, py::arg("n"),
+             "Select the n minimal elements.")
         //
         .def("__len__", [](ten_ptr self) { return self->shape()[0]; })
         .def("__repr__", &Tensor___repr__)
@@ -508,38 +531,40 @@ void export_Tensor(py::module& m) {
         //
         .def("__pos__", [](ten_ptr self) { return self; })               // + tensor
         .def("__neg__", [](ten_ptr self) { return self->scale(-1.0); })  // - tensor
-        .def("__add__", &Tensor_scalar__add__)            // tensor + scalar
-        .def("__sub__", &Tensor_scalar__sub__)            // tensor - scalar
-        .def("__radd__", &Tensor_scalar__add__)           // scalar + tensor
-        .def("__rsub__", &Tensor_scalar__sub__)           // scalar - tensor
-        .def("__imul__", &Tensor_scalar__imul__)          // tensor *= scalar
-        .def("__mul__", &Tensor_scalar__mul__)            // tensor * scalar
-        .def("__rmul__", &Tensor_scalar__mul__)           // scalar * tensor
-        .def("__itruediv__", &Tensor_scalar__itruediv__)  // tensor /= scalar
-        .def("__truediv__", &Tensor_scalar__truediv__)    // tensor / scalar
-                                                          //
-        .def("__mul__", &Tensor__mul__,
+        .def("__add__", &Tensor_scalar__add__, py::arg("number"))    // tensor + scalar
+        .def("__sub__", &Tensor_scalar__sub__, py::arg("number"))    // tensor - scalar
+        .def("__radd__", &Tensor_scalar__add__, py::arg("number"))   // scalar + tensor
+        .def("__rsub__", &Tensor_scalar__sub__, py::arg("number"))   // scalar - tensor
+        .def("__imul__", &Tensor_scalar__imul__, py::arg("number"))  // tensor *= scalar
+        .def("__mul__", &Tensor_scalar__mul__, py::arg("number"))    // tensor * scalar
+        .def("__rmul__", &Tensor_scalar__mul__, py::arg("number"))   // scalar * tensor
+        .def("__itruediv__", &Tensor_scalar__itruediv__,
+             py::arg("number"))  // tensor /= scalar
+        .def("__truediv__", &Tensor_scalar__truediv__,
+             py::arg("number"))  // tensor / scalar
+                                 //
+        .def("__mul__", &Tensor__mul__, py::arg("other"),
              "Multiply two tensors elementwise.")  // tensor * tensor
-        .def("__truediv__", &Tensor__truediv__,
-             "Divide two tensors elementwise.")  // tensor / tensor
-        .def("__iadd__", &Tensor__iadd__)        // tensor += tensor
-        .def("__add__", &Tensor__add__)          // tensor + tensor
-        .def("__isub__", &Tensor__isub__)        // tensor -= tensor
-        .def("__sub__", &Tensor__sub__)          // tensor - tensor
+        .def("__truediv__", &Tensor__truediv__, py::arg("other"),
+             "Divide two tensors elementwise.")              // tensor / tensor
+        .def("__iadd__", &Tensor__iadd__, py::arg("other"))  // tensor += tensor
+        .def("__add__", &Tensor__add__, py::arg("other"))    // tensor + tensor
+        .def("__isub__", &Tensor__isub__, py::arg("other"))  // tensor -= tensor
+        .def("__sub__", &Tensor__sub__, py::arg("other"))    // tensor - tensor
         //
-        .def("__matmul__", &Tensor__matmul__)  // tensor @ tensor
+        .def("__matmul__", &Tensor__matmul__, py::arg("other"))  // tensor @ tensor
         //
         ;
 
-  m.def("evaluate", &evaluate);
-  m.def("tensordot", &tensordot_1, "a"_a, "b"_a, "axes"_a);
-  m.def("tensordot", &tensordot_2, "a"_a, "b"_a, "axes"_a);
-  m.def("tensordot", &tensordot_3, "a"_a, "b"_a);
-  m.def("direct_sum", &direct_sum, "a"_a, "b"_a);
-  m.def("trace", &Tensor_trace_1, "subscripts"_a, "tensor"_a);
-  m.def("trace", &Tensor_trace_2, "tensor"_a);
-  m.def("linear_combination_strict", &linear_combination_strict, "coefficients"_a,
-        "tensors"_a);
+  m.def("evaluate", &evaluate, py::arg("tensor"));
+  m.def("tensordot", &tensordot_1, py::arg("a"), py::arg("b"), py::arg("axes"));
+  m.def("tensordot", &tensordot_2, py::arg("a"), py::arg("b"), py::arg("axes"));
+  m.def("tensordot", &tensordot_3, py::arg("a"), py::arg("b"));
+  m.def("direct_sum", &direct_sum, py::arg("a"), py::arg("b"));
+  m.def("trace", &Tensor_trace_1, py::arg("subscripts"), py::arg("tensor"));
+  m.def("trace", &Tensor_trace_2, py::arg("tensor"));
+  m.def("linear_combination_strict", &linear_combination_strict, py::arg("coefficients"),
+        py::arg("tensors"));
 }
 
 }  // namespace libadcc
diff --git a/libadcc_src/pyiface/export_adc_pp.cc b/libadcc_src/pyiface/export_adc_pp.cc
index 784da441..a5990ce3 100644
--- a/libadcc_src/pyiface/export_adc_pp.cc
+++ b/libadcc_src/pyiface/export_adc_pp.cc
@@ -25,21 +25,22 @@
 namespace libadcc {
 
 namespace py = pybind11;
-using namespace pybind11::literals;
 
 void export_adc_pp(py::module& m) {
   m.def("amplitude_vector_enforce_spin_kind", &amplitude_vector_enforce_spin_kind,
+        py::arg("doubles_tensor"), py::arg("block"), py::arg("spin_kind"),
         "Apply the spin symmetrisation required to make the doubles and higher parts of "
         "an amplitude vector consist of components for a particular spin kind only.");
 
-  m.def("fill_pp_doubles_guesses", &fill_pp_doubles_guesses, "guesses_d"_a, "mospaces"_a,
-        "df02"_a, "df13"_a, "spin_change_twice"_a, "degeneracy_tolerance"_a,
+  m.def("fill_pp_doubles_guesses", &fill_pp_doubles_guesses, py::arg("guesses_d"),
+        py::arg("mospaces"), py::arg("df1"), py::arg("df2"), py::arg("spin_change_twice"),
+        py::arg("degeneracy_tolerance"),
         "Fill the passed vector of doubles blocks with doubles guesses using the "
-        "delta-Fock matrices df02 and df13, which are the two delta-Fock matrices "
+        "delta-Fock matrices df1 and df2, which are the two delta-Fock matrices "
         "involved in the doubles block.\n\nguesses_d    Vectors of guesses, all elements "
         "are assumed to be initialised to zero and the symmetry is assumed to be "
-        "properly set up.\nmospaces     Mospaces object\ndf02         Delta-Fock between "
-        "spaces 0 and 2 of the ADC matrix\ndf13         Delta-Fock between spaces 1 and "
+        "properly set up.\nmospaces     Mospaces object\ndf1         Delta-Fock between "
+        "spaces 0 and 2 of the ADC matrix\ndf2         Delta-Fock between spaces 1 and "
         "3 of the ADC matrix\nspin_change_twice   Twice the value of the spin change to "
         "enforce in an excitation.\ndegeneracy_tolerance  Tolerance for two entries of "
         "the diagonal to be considered degenerate, i.e. identical.\nReturns     The "
diff --git a/libadcc_src/pyiface/export_threading.cc b/libadcc_src/pyiface/export_threading.cc
index 996e584c..ccbaeeaa 100644
--- a/libadcc_src/pyiface/export_threading.cc
+++ b/libadcc_src/pyiface/export_threading.cc
@@ -40,7 +40,7 @@ void export_threading(py::module& m) {
   m.def(
         "get_n_threads", [threadpool_ptr]() { return threadpool_ptr->n_running(); },
         "Get the number of running worker threads used by adcc.");
-  m.def("set_n_threads", set_threads,
+  m.def("set_n_threads", set_threads, py::arg("n_threads"),
         "Set the number of running worker threads used by adcc");
   m.def(
         "set_n_threads_total",
@@ -48,6 +48,7 @@ void export_threading(py::module& m) {
           const size_t n_running = threadpool_ptr->n_running();
           threadpool_ptr->reinit(n_running, n_total);
         },
+        py::arg("n_total"),
         "Set the total number of threads (running and sleeping) used by adcc. This will "
         "disappear in the future. Do not rely on it.");
   m.def(
diff --git a/libadcc_src/pyiface/ndarray.hh b/libadcc_src/pyiface/ndarray.hh
new file mode 100644
index 00000000..4a10ee47
--- /dev/null
+++ b/libadcc_src/pyiface/ndarray.hh
@@ -0,0 +1,96 @@
+//
+// Copyright (C) 2026 by the adcc authors
+//
+// This file is part of adcc.
+//
+// adcc is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// adcc is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with adcc. If not, see <http://www.gnu.org/licenses/>.
+//
+
+#pragma once
+#include "../exceptions.hh"
+#include <cstddef>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <string>
+
+namespace libadcc {
+
+namespace py = pybind11;
+
+// Small wrapper around py::array_t to enable the correct
+// type hints for the dimensionality of numpy arrays in the stub file.
+// The dtype is handled exactly like in py::array_t<T>, i.e. other types are
+// silently casted (and copied) to T if necessary (e.g. int -> double).
+// The dimensionality on the other hand is never adapted: constructing from an
+// array with a different ndim throws.
+template <typename T, size_t Ndim>
+class NDArray : public py::array_t<T> {
+  static_assert(Ndim > 0, "NDArray requires at least one dimension");
+
+ public:
+  // Constructors:
+  // Needed by the type caster: it default-constructs via reinterpret_steal and
+  // loads via reinterpret_borrow (after check_ passed), so both are required as
+  // soon as NDArray appears as a function *parameter* (fill_*).
+  NDArray(py::handle h, py::object::borrowed_t b) : py::array_t<T>(h, b) {}
+  NDArray(py::handle h, py::object::stolen_t s) : py::array_t<T>(h, s) {}
+
+  // Converting constructor.
+  NDArray(py::array array) : py::array_t<T>(validate(array)) {}
+
+  // Hides py::array_t<T>::check_, which only validates the dtype. isinstance<NDArray>
+  // dispatches to this, i.e., the type caster (which builds the value through
+  // reinterpret_borrow and therefore bypasses the converting constructor above)
+  // rejects arrays of the wrong dimensionality.
+  static bool check_(py::handle h) {
+    return py::array_t<T>::check_(h) &&
+           py::reinterpret_borrow<py::array>(h).ndim() == static_cast<py::ssize_t>(Ndim);
+  }
+
+ private:
+  static const py::array& validate(const py::array& array) {
+    if (array.ndim() != static_cast<py::ssize_t>(Ndim)) {
+      throw dimension_mismatch("Expected a " + std::to_string(Ndim) +
+                               "-dimensional array, but got " +
+                               std::to_string(array.ndim()) + " dimensions.");
+    }
+    return array;
+  }
+};
+}  // namespace libadcc
+
+namespace pybind11 {
+namespace detail {
+
+// Create a string of Ndim "int, int, ..." as shape for the np.ndarray type hint
+template <size_t Ndim>
+constexpr auto ndim_name() {
+  // ndim_name<0> would underflow to ndim_name<SIZE_MAX> (Ndim == 1 is specialized below)
+  static_assert(Ndim > 0, "NDArray requires at least one dimension");
+  return ndim_name<Ndim - 1>() + const_name(", int");
+}
+template <>
+constexpr auto ndim_name<1>() {
+  return const_name("int");
+}
+
+template <typename T, std::size_t Ndim>
+struct handle_type_name<libadcc::NDArray<T, Ndim>> {
+  static constexpr auto name = const_name("numpy.ndarray[tuple[") + ndim_name<Ndim>() +
+                               const_name("], numpy.dtype[") +
+                               npy_format_descriptor<T>::name + const_name("]]");
+};
+
+}  // namespace detail
+}  // namespace pybind11
diff --git a/libadcc_src/pyiface/util.cc b/libadcc_src/pyiface/util.cc
index f6f750e8..d6d60f26 100644
--- a/libadcc_src/pyiface/util.cc
+++ b/libadcc_src/pyiface/util.cc
@@ -18,16 +18,14 @@
 //
 
 #include "util.hh"
-#include "../config.hh"
 #include "../exceptions.hh"
-#include <sstream>
-#include <type_traits>
+#include "pybind11/typing.h"
 
 namespace libadcc {
 
 namespace py = pybind11;
 
-py::tuple shape_tuple(const std::vector<size_t>& shape) {
+py::typing::Tuple<size_t, py::ellipsis> shape_tuple(const std::vector<size_t>& shape) {
   switch (shape.size()) {
     case 0:
       throw runtime_error("Encountered unexpected dimensionality 0.");
diff --git a/libadcc_src/pyiface/util.hh b/libadcc_src/pyiface/util.hh
index 26573d59..aab1f356 100644
--- a/libadcc_src/pyiface/util.hh
+++ b/libadcc_src/pyiface/util.hh
@@ -22,6 +22,7 @@
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
+#include <pybind11/typing.h>
 #include <vector>
 
 namespace libadcc {
@@ -29,7 +30,7 @@ namespace libadcc {
 namespace py = pybind11;
 
 /** Make a py::tuple from a vector representing the shape */
-py::tuple shape_tuple(const std::vector<size_t>& shape);
+py::typing::Tuple<size_t, py::ellipsis> shape_tuple(const std::vector<size_t>& shape);
 
 /** Convert a list of tensors to a vector of shared pointers to Tensor */
 template <typename Listlike>
diff --git a/libadcc_src/tests/HFSolutionMock.hh b/libadcc_src/tests/HFSolutionMock.hh
index 5b2f20f6..d7fd4a4f 100644
--- a/libadcc_src/tests/HFSolutionMock.hh
+++ b/libadcc_src/tests/HFSolutionMock.hh
@@ -49,11 +49,13 @@ struct HFSolutionMock : public HartreeFockSolution_i {
   size_t spin_multiplicity() const override { return restricted() ? 1 : 0; }
 
   size_t n_bas() const override { return exposed_n_bas; }
-  void nuclear_multipole(size_t /*order*/, std::array<scalar_type, 3> /*gauge_origin*/,
-                         scalar_type* /*buffer*/, size_t /*size*/) const override {
+  void nuclear_multipole(
+        size_t /*order*/,
+        std::tuple<scalar_type, scalar_type, scalar_type> /*gauge_origin*/,
+        scalar_type* /*buffer*/, size_t /*size*/) const override {
     throw not_implemented_error("Not implemented.");
   }
-  const std::array<scalar_type, 3> gauge_origin_to_xyz(
+  const std::tuple<scalar_type, scalar_type, scalar_type> gauge_origin_to_xyz(
         std::string /*gauge_origin*/) const override {
     throw not_implemented_error("Not implemented.");
   }
diff --git a/pyproject.toml b/pyproject.toml
index feb37925..e16c45cb 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -47,7 +47,7 @@ tests = ["pytest", "pytest-cov", "pandas >= 0.25.0"]
 # since we explicity call setup.py
 cpptest = [
 	"setuptools >= 77.0",
-	"pybind11 >= 2.6"
+	"pybind11 >= 3.0"
 ]
 build_docs = [
 	"sphinx>=2", "breathe", "sphinxcontrib-bibtex",
@@ -64,7 +64,7 @@ Homepage = "https://adc-connect.org"
 [build-system]
 requires = [
 	"setuptools >= 77.0",  # due to license
-	"pybind11 >= 2.6"
+	"pybind11 >= 3.0"
 ]
 build-backend = "setuptools.build_meta"
 
diff --git a/scripts/generate_libadcc_stub.sh b/scripts/generate_libadcc_stub.sh
index dc2878b0..0c9f2b56 100755
--- a/scripts/generate_libadcc_stub.sh
+++ b/scripts/generate_libadcc_stub.sh
@@ -4,11 +4,14 @@ set -e
 cd ..
 pybind11-stubgen libadcc -o .
 # Remove the __backend__ module-level variable added by ExportAdcc.cc
+# Unions currently produce 'libadcc.Tensor | T', which is not correctly
+# resolved to 'Tensor | T' by pybind11-stubgen.
 python3 -c """
 import re, pathlib
 p = pathlib.Path('libadcc.pyi')
 text = p.read_text()
 text = re.sub(r'\n__backend__\s*:.*?}\n', '', text, flags=re.DOTALL)
+text = text.replace('libadcc.Tensor', 'Tensor')
 p.write_text(text)
 """
 ruff format libadcc.pyi