diff --git a/CMakeLists.txt b/CMakeLists.txt
index 33d7b0dcd5..d53eb20ed3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -332,6 +332,7 @@ set(SeQuant_src
         SeQuant/core/runtime.hpp
         SeQuant/core/space.cpp
         SeQuant/core/space.hpp
+        SeQuant/core/slotted_index.hpp
         SeQuant/core/tag.hpp
         SeQuant/core/tensor_canonicalizer.cpp
         SeQuant/core/tensor_canonicalizer.hpp
@@ -359,6 +360,7 @@ set(SeQuant_src
         SeQuant/core/utility/macros.hpp
         SeQuant/core/utility/memoize.hpp
         SeQuant/core/utility/nodiscard.hpp
+        SeQuant/core/utility/overloads.hpp
         SeQuant/core/utility/permutation.hpp
         SeQuant/core/utility/scope.hpp
         SeQuant/core/utility/singleton.hpp
@@ -381,6 +383,7 @@ set(SeQuant_src
         SeQuant/domain/mbpt/context.hpp
         SeQuant/domain/mbpt/convention.cpp
         SeQuant/domain/mbpt/convention.hpp
+        SeQuant/domain/mbpt/detail/concepts.hpp
         SeQuant/domain/mbpt/models/cc.cpp
         SeQuant/domain/mbpt/models/cc.hpp
         SeQuant/domain/mbpt/op_registry.hpp
diff --git a/SeQuant/core/attr.hpp b/SeQuant/core/attr.hpp
index a83e3d4067..5d8e499a17 100644
--- a/SeQuant/core/attr.hpp
+++ b/SeQuant/core/attr.hpp
@@ -9,6 +9,7 @@
 
 #include <cassert>
 #include <cstdlib>
+#include <ostream>
 #include <string>
 
 namespace sequant {
@@ -79,6 +80,36 @@ inline std::wstring to_wolfram(BraKetPos a) {
   }
 }
 
+/// index slot types
+///
+/// @note This does not include slot bundles, like braket, etc.
+enum class SlotType {
+  Bra,
+  Ket,
+  Aux,
+  Proto,
+};
+
+template <typename CharT, typename Traits>
+std::basic_ostream<CharT, Traits>& operator<<(
+    std::basic_ostream<CharT, Traits>& stream, SlotType origin) {
+  switch (origin) {
+    case SlotType::Bra:
+      stream << "Bra";
+      break;
+    case SlotType::Ket:
+      stream << "Ket";
+      break;
+    case SlotType::Aux:
+      stream << "Aux";
+      break;
+    case SlotType::Proto:
+      stream << "Proto";
+      break;
+  }
+  return stream;
+}
+
 enum class Statistics {
   FermiDirac,
   BoseEinstein,
diff --git a/SeQuant/core/expressions/abstract_tensor.hpp b/SeQuant/core/expressions/abstract_tensor.hpp
index a41699d465..329381fd1a 100644
--- a/SeQuant/core/expressions/abstract_tensor.hpp
+++ b/SeQuant/core/expressions/abstract_tensor.hpp
@@ -15,7 +15,6 @@
 
 #include <cstdlib>
 #include <memory>
-#include <ostream>
 #include <stdexcept>
 #include <string>
 #include <string_view>
@@ -27,31 +26,6 @@
 
 namespace sequant {
 
-/// index slot types
-///
-/// @note This does not include slot bundles, like braket, etc.
-enum class SlotType { Bra = 0, Ket = 1, Aux = 2, Proto = 3 };
-
-template <typename CharT, typename Traits>
-std::basic_ostream<CharT, Traits>& operator<<(
-    std::basic_ostream<CharT, Traits>& stream, SlotType origin) {
-  switch (origin) {
-    case SlotType::Bra:
-      stream << "Bra";
-      break;
-    case SlotType::Ket:
-      stream << "Ket";
-      break;
-    case SlotType::Aux:
-      stream << "Aux";
-      break;
-    case SlotType::Proto:
-      stream << "Proto";
-      break;
-  }
-  return stream;
-}
-
 class TensorCanonicalizer;
 
 /// AbstractTensor is a [tensor](https://en.wikipedia.org/wiki/Tensor) over
diff --git a/SeQuant/core/expressions/result_expr.hpp b/SeQuant/core/expressions/result_expr.hpp
index bc26db069d..2e67cfb5b3 100644
--- a/SeQuant/core/expressions/result_expr.hpp
+++ b/SeQuant/core/expressions/result_expr.hpp
@@ -6,6 +6,7 @@
 #include <SeQuant/core/expressions/tensor.hpp>
 #include <SeQuant/core/expressions/variable.hpp>
 #include <SeQuant/core/index.hpp>
+#include <SeQuant/core/slotted_index.hpp>
 #include <SeQuant/core/utility/macros.hpp>
 
 #include <initializer_list>
@@ -102,15 +103,24 @@ class ResultExpr {
     // based on the particle they belong to and that bra and
     // ket indices are assigned to the same set of particles.
     for (std::size_t i = 0; i < m_braIndices.size(); ++i) {
-      if constexpr (std::is_constructible_v<Group,
-                                            std::initializer_list<Index>>) {
+      if constexpr (std::is_constructible_v<
+                        Group, std::initializer_list<SlottedIndex>>) {
+        groups.emplace_back(std::initializer_list<SlottedIndex>{
+            {m_braIndices.at(i), SlotType::Bra},
+            {m_ketIndices.at(i), SlotType::Ket}});
+      } else if constexpr (std::is_constructible_v<
+                               Group, std::initializer_list<Index>>) {
         groups.emplace_back(std::initializer_list<Index>{m_braIndices.at(i),
                                                          m_ketIndices.at(i)});
+      } else if constexpr (std::is_constructible_v<Group, SlottedIndex,
+                                                   SlottedIndex>) {
+        groups.emplace_back(SlottedIndex{m_braIndices.at(i), SlotType::Bra},
+                            SlottedIndex{m_ketIndices.at(i), SlotType::Ket});
       } else {
         static_assert(
             std::is_constructible_v<Group, Index, Index>,
-            "Group is expected to be constructible from two indices or from an "
-            "initializer_list of indices");
+            "Group is expected to be constructible from two (slotted) indices "
+            "or from an initializer_list of (slotted) indices");
         groups.emplace_back(m_braIndices.at(i), m_ketIndices.at(i));
       }
     }
diff --git a/SeQuant/core/meta.hpp b/SeQuant/core/meta.hpp
index 23b20f76bb..8de5139f24 100644
--- a/SeQuant/core/meta.hpp
+++ b/SeQuant/core/meta.hpp
@@ -7,9 +7,10 @@
 
 #include <complex>
 #include <memory>
+#include <type_traits>
+
 #include <range/v3/range/access.hpp>
 #include <range/v3/range/traits.hpp>
-#include <type_traits>
 
 namespace sequant {
 
@@ -506,6 +507,81 @@ struct std_array_size<std::array<T, N>>
 template <typename T>
 constexpr inline std::size_t std_array_size_v = std_array_size<T>::value;
 
+/// make_immutable_t
+/// Contrary to std::add_const_t, this works as one would expect on reference
+/// types
+template <typename T>
+using make_immutable_t = std::conditional_t<
+    std::is_lvalue_reference_v<T>,
+    std::add_lvalue_reference_t<std::add_const_t<std::remove_reference_t<T>>>,
+    std::conditional_t<std::is_rvalue_reference_v<T>,
+                       std::add_rvalue_reference_t<
+                           std::add_const_t<std::remove_reference_t<T>>>,
+                       std::add_const_t<T>>>;
+static_assert(std::is_const_v<make_immutable_t<float>>);
+static_assert(std::is_const_v<make_immutable_t<const float>>);
+static_assert(std::is_lvalue_reference_v<make_immutable_t<float &>>);
+static_assert(std::is_lvalue_reference_v<make_immutable_t<const float &>>);
+static_assert(
+    std::is_const_v<std::remove_reference_t<make_immutable_t<float &>>>);
+static_assert(
+    std::is_const_v<std::remove_reference_t<make_immutable_t<const float &>>>);
+static_assert(std::is_rvalue_reference_v<make_immutable_t<float &&>>);
+static_assert(std::is_rvalue_reference_v<make_immutable_t<const float &&>>);
+static_assert(
+    std::is_const_v<std::remove_reference_t<make_immutable_t<float &&>>>);
+static_assert(
+    std::is_const_v<std::remove_reference_t<make_immutable_t<const float &&>>>);
+
+/// make_mutable_t
+/// Contrary to std::remove_const_t, this works as one would expect on reference
+/// types
+template <typename T>
+using make_mutable_t = std::conditional_t<
+    std::is_lvalue_reference_v<T>,
+    std::add_lvalue_reference_t<
+        std::remove_const_t<std::remove_reference_t<T>>>,
+    std::conditional_t<std::is_rvalue_reference_v<T>,
+                       std::add_rvalue_reference_t<
+                           std::remove_const_t<std::remove_reference_t<T>>>,
+                       std::remove_const_t<T>>>;
+static_assert(!std::is_const_v<make_mutable_t<float>>);
+static_assert(!std::is_const_v<make_mutable_t<const float>>);
+static_assert(std::is_lvalue_reference_v<make_mutable_t<float &>>);
+static_assert(std::is_lvalue_reference_v<make_mutable_t<const float &>>);
+static_assert(
+    !std::is_const_v<std::remove_reference_t<make_mutable_t<float &>>>);
+static_assert(
+    !std::is_const_v<std::remove_reference_t<make_mutable_t<const float &>>>);
+static_assert(std::is_rvalue_reference_v<make_mutable_t<float &&>>);
+static_assert(std::is_rvalue_reference_v<make_mutable_t<const float &&>>);
+static_assert(
+    !std::is_const_v<std::remove_reference_t<make_mutable_t<float &&>>>);
+static_assert(
+    !std::is_const_v<std::remove_reference_t<make_mutable_t<const float &&>>>);
+
+/// is_immutable_v
+/// Contrary to std::is_const_v, this works as one would expect on reference
+/// types
+template <typename T>
+constexpr bool is_immutable_v = std::is_const_v<std::remove_reference_t<T>>;
+static_assert(is_immutable_v<const float>);
+static_assert(!is_immutable_v<float>);
+static_assert(is_immutable_v<const float &>);
+static_assert(!is_immutable_v<float &>);
+static_assert(is_immutable_v<const float &&>);
+static_assert(!is_immutable_v<float &&>);
+
+/// mimic_constness_t
+/// Makes To const, if From is const, else makes To non-const
+template <typename From, typename To>
+using mimic_constness_t =
+    std::conditional_t<is_immutable_v<From>, make_immutable_t<To>,
+                       make_mutable_t<To>>;
+static_assert(
+    std::same_as<mimic_constness_t<const int &, float &>, const float &>);
+static_assert(std::same_as<mimic_constness_t<int &, const float &>, float &>);
+
 }  // namespace meta
 }  // namespace sequant
 
diff --git a/SeQuant/core/slotted_index.hpp b/SeQuant/core/slotted_index.hpp
new file mode 100644
index 0000000000..4c655813df
--- /dev/null
+++ b/SeQuant/core/slotted_index.hpp
@@ -0,0 +1,42 @@
+#ifndef SEQUANT_SLOTTED_INDEX_H
+#define SEQUANT_SLOTTED_INDEX_H
+
+#include <SeQuant/core/attr.hpp>
+#include <SeQuant/core/index.hpp>
+
+namespace sequant {
+
+/// Combination of an Index along with information about the Slot it occupies
+///
+/// @note Usage of this class only makes sense, if the index in question does
+/// not appear in multiple slots
+///
+/// @sa Index
+/// @sa Slot
+class SlottedIndex {
+ public:
+  SlottedIndex(Index idx, SlotType slot)
+      : idx_(std::move(idx)), slot_(std::move(slot)) {}
+
+  Index &index() { return idx_; }
+
+  const Index &index() const { return idx_; }
+
+  SlotType slot_type() const { return slot_; }
+
+  friend bool operator==(const SlottedIndex &lhs, const SlottedIndex &rhs) {
+    return lhs.index() == rhs.index() && lhs.slot_type() == rhs.slot_type();
+  }
+
+  friend bool operator!=(const SlottedIndex &lhs, const SlottedIndex &rhs) {
+    return !(lhs == rhs);
+  }
+
+ private:
+  Index idx_;
+  SlotType slot_;
+};
+
+}  // namespace sequant
+
+#endif
diff --git a/SeQuant/core/utility/indices.hpp b/SeQuant/core/utility/indices.hpp
index 504ca00661..e348547ac1 100644
--- a/SeQuant/core/utility/indices.hpp
+++ b/SeQuant/core/utility/indices.hpp
@@ -1,20 +1,25 @@
 #ifndef SEQUANT_CORE_UTILITY_INDICES_HPP
 #define SEQUANT_CORE_UTILITY_INDICES_HPP
 
+#include <SeQuant/core/attr.hpp>
 #include <SeQuant/core/container.hpp>
 #include <SeQuant/core/expr.hpp>
 #include <SeQuant/core/index.hpp>
+#include <SeQuant/core/meta.hpp>
 #include <SeQuant/core/op.hpp>
 #include <SeQuant/core/reserved.hpp>
+#include <SeQuant/core/slotted_index.hpp>
 #include <SeQuant/core/utility/macros.hpp>
 
 #include <range/v3/view.hpp>
 
 #include <algorithm>
+#include <concepts>
 #include <iterator>
 #include <optional>
 #include <ranges>
 #include <set>
+#include <tuple>
 #include <type_traits>
 #include <vector>
 
@@ -434,8 +439,11 @@ std::string csv_labels(Rng&& idxs) {
 ///
 /// @see get_bra_idx
 /// @see get_ket_idx
-template <typename Container = container::svector<container::svector<Index>>>
-Container external_indices(const Expr& expr) {
+template <template <class> class Container = container::svector,
+          template <class> class Group = container::svector>
+Container<Group<SlottedIndex>> external_indices(const Expr& expr) {
+  using HolderType = Container<Group<SlottedIndex>>;
+
   if (!expr.is<Sum>() && !expr.is<Product>() && !expr.is<Tensor>()) {
     return {};
   }
@@ -445,16 +453,17 @@ Container external_indices(const Expr& expr) {
 
     const std::size_t num_braket =
         std::max(tensor.bra_rank(), tensor.ket_rank());
-    Container cont(num_braket + tensor.aux_rank());
+    HolderType cont;
+    cont.resize(num_braket + tensor.aux_rank());
 
     for (std::size_t i = 0; i < tensor.bra_rank(); ++i) {
-      cont.at(i).push_back(tensor.bra()[i]);
+      cont.at(i).emplace_back(tensor.bra()[i], SlotType::Bra);
     }
     for (std::size_t i = 0; i < tensor.ket_rank(); ++i) {
-      cont.at(i).push_back(tensor.ket()[i]);
+      cont.at(i).emplace_back(tensor.ket()[i], SlotType::Ket);
     }
     for (std::size_t i = 0; i < tensor.aux_rank(); ++i) {
-      cont.at(i + num_braket).push_back(tensor.aux()[i]);
+      cont.at(i + num_braket).emplace_back(tensor.aux()[i], SlotType::Aux);
     }
 
     if (tensor.label() == reserved::symm_label() ||
@@ -463,7 +472,7 @@ Container external_indices(const Expr& expr) {
       // indices are conjugated (reversed). Hence, we have to swap the
       // determined bra and ket indices.
       for (std::size_t i = 0; i < num_braket; ++i) {
-        std::swap(cont.at(i).at(0), cont.at(i).at(1));
+        std::swap(cont.at(i).at(0).index(), cont.at(i).at(1).index());
       }
     }
 
@@ -485,26 +494,24 @@ Container external_indices(const Expr& expr) {
 
   if (symmetrizer.has_value()) {
     // Generate external index list from symmetrization operator
-    // However, the symmetrizer has bra/ket conjugated (reversed)
-    Container res = external_indices<Container>(symmetrizer.value());
-
-    return res;
+    return external_indices<Container, Group>(symmetrizer.value());
   }
 
   IndexGroups groups = get_unique_indices<container::svector<Index>>(expr);
 
   const std::size_t num_braket = std::max(groups.bra.size(), groups.ket.size());
 
-  Container cont(num_braket + groups.aux.size());
+  HolderType cont;
+  cont.resize(num_braket + groups.aux.size());
 
   for (std::size_t i = 0; i < groups.bra.size(); ++i) {
-    cont.at(i).push_back(groups.bra[i]);
+    cont.at(i).emplace_back(groups.bra[i], SlotType::Bra);
   }
   for (std::size_t i = 0; i < groups.ket.size(); ++i) {
-    cont.at(i).push_back(groups.ket[i]);
+    cont.at(i).emplace_back(groups.ket[i], SlotType::Ket);
   }
   for (std::size_t i = 0; i < groups.aux.size(); ++i) {
-    cont.at(num_braket + i).push_back(groups.aux[i]);
+    cont.at(num_braket + i).emplace_back(groups.aux[i], SlotType::Aux);
   }
 
   return cont;
@@ -521,55 +528,185 @@ Container external_indices(const Expr& expr) {
 ///
 /// @see get_bra_idx
 /// @see get_ket_idx
-template <typename Container = container::svector<container::svector<Index>>>
-Container external_indices(const ExprPtr& expr) {
-  return external_indices<Container>(*expr);
+template <template <class> class Container = container::svector,
+          template <class> class Group = container::svector>
+Container<Group<SlottedIndex>> external_indices(const ExprPtr& expr) {
+  return external_indices<Container, Group>(*expr);
 }
 
-/// @param container An index group expected to represent a pair of indices
-/// @returns The bra index of the provided pair
-///
-/// @tparam T The type of the container used to represent the index pair. May
-/// be a std::pair<Index, Index> or a container that retains the order of its
-/// elements based on how they have been inserted.
 template <typename T>
-auto get_bra_idx(T&& container)
-    -> std::conditional_t<std::is_const_v<std::remove_reference_t<T>>,
-                          const Index&, Index&> {
-  if constexpr (std::is_same_v<std::remove_cvref_t<T>,
-                               std::pair<Index, Index>>) {
-    return container.first;
+concept SlottedIndexContainer =
+    std::ranges::range<T> &&
+    std::same_as<std::ranges::range_value_t<T>, SlottedIndex>;
+
+template <typename T>
+concept IndexContainer =
+    std::ranges::range<T> && std::same_as<std::ranges::range_value_t<T>, Index>;
+
+template <typename T>
+concept SlottedIndexTuple =
+    std::is_same_v<std::tuple_element_t<0, std::remove_cvref_t<T>>,
+                   SlottedIndex> &&
+    std::is_same_v<std::tuple_element_t<1, std::remove_cvref_t<T>>,
+                   SlottedIndex>;
+template <typename T>
+concept IndexTuple =
+    std::is_same_v<std::tuple_element_t<0, std::remove_cvref_t<T>>, Index> &&
+    std::is_same_v<std::tuple_element_t<1, std::remove_cvref_t<T>>, Index>;
+
+template <typename T>
+concept SlottedIndexGroup = SlottedIndexContainer<T> || SlottedIndexTuple<T>;
+
+template <typename T>
+concept IndexGroup = IndexContainer<T> || IndexTuple<T>;
+
+template <typename T>
+concept SlottedIndexGroupContainer =
+    std::ranges::range<T> && SlottedIndexGroup<std::ranges::range_value_t<T>>;
+
+template <typename T>
+concept IndexGroupContainer =
+    std::ranges::range<T> && IndexGroup<std::ranges::range_value_t<T>>;
+
+/// @param group A group of indices with slot information
+/// @returns The bra index of the provided group
+///
+///@note This function assumes that there are no duplicate slots
+///
+/// @tparam Group The type of the container used to represent the index group.
+/// May be a std::pair<SlottedIndex, SlottedIndex> or a container that retains
+/// the order of its elements based on how they have been inserted
+/// (vector-like).
+template <SlottedIndexGroup Group>
+decltype(auto) get_bra_idx(Group&& group) {
+  if constexpr (SlottedIndexTuple<Group>) {
+    SEQUANT_ASSERT(std::get<0>(group).slot_type() == SlotType::Bra ||
+                   std::get<1>(group).slot_type() == SlotType::Bra);
+    return std::get<0>(group).slot_type() == SlotType::Bra
+               ? std::get<0>(group).index()
+               : std::get<1>(group).index();
   } else {
     using std::begin;
-    using std::size;
-    SEQUANT_ASSERT(size(container) == 2);
-    return *begin(container);
+    using std::end;
+    auto it =
+        std::find_if(begin(group), end(group), [](const SlottedIndex& idx) {
+          return idx.slot_type() == SlotType::Bra;
+        });
+    SEQUANT_ASSERT(it != end(group));
+    return it->index();
   }
 }
 
-/// @param container An index group expected to represent a pair of indices
-/// @returns The ket index of the provided pair
+/// @param group A group of indices assumed to be in canonical ordering
+/// @returns The bra index of the provided group
 ///
-/// @tparam T The type of the container used to represent the index pair. May
-/// be a std::pair<Index, Index> or a container that retains the order of its
-/// elements based on how they have been inserted.
-template <typename T>
-auto get_ket_idx(T&& container)
-    -> std::conditional_t<std::is_const_v<std::remove_reference_t<T>>,
-                          const Index&, Index&> {
-  if constexpr (std::is_same_v<std::remove_cvref_t<T>,
-                               std::pair<Index, Index>>) {
-    return container.second;
+/// @tparam Group The type of the container used to represent the index group.
+/// May be a std::pair<Index, Index> or a container that retains
+/// the order of its elements based on how they have been inserted
+/// (vector-like).
+template <IndexGroup Group>
+decltype(auto) get_bra_idx(Group&& group) {
+  if constexpr (IndexTuple<Group>) {
+    return std::get<0>(group);
   } else {
     using std::begin;
-    using std::size;
-    SEQUANT_ASSERT(size(container) == 2);
-    auto it = begin(container);
+    using std::end;
+    auto it = begin(group);
+    SEQUANT_ASSERT(it != end(group));
+    return *it;
+  }
+}
+
+/// @param group A group of indices with slot information
+/// @returns The ket index of the provided group
+///
+///@note This function assumes that there are no duplicate slots
+///
+/// @tparam Group The type of the container used to represent the index group.
+/// May be a std::pair<SlottedIndex, SlottedIndex> or a container that retains
+/// the order of its elements based on how they have been inserted
+/// (vector-like).
+template <SlottedIndexGroup Group>
+decltype(auto) get_ket_idx(Group&& group) {
+  if constexpr (SlottedIndexTuple<Group>) {
+    SEQUANT_ASSERT(std::get<0>(group).slot_type() == SlotType::Ket ||
+                   std::get<1>(group).slot_type() == SlotType::Ket);
+    return std::get<1>(group).slot_type() == SlotType::Ket
+               ? std::get<1>(group).index()
+               : std::get<0>(group).index();
+  } else {
+    // Note: We're using reverse iteration order as typically, we expect the ket
+    // slot to be the second of two
+    using std::rbegin;
+    using std::rend;
+    auto it =
+        std::find_if(rbegin(group), rend(group), [](const SlottedIndex& idx) {
+          return idx.slot_type() == SlotType::Ket;
+        });
+    SEQUANT_ASSERT(it != rend(group));
+    return it->index();
+  }
+}
+
+/// @param group A group of indices assumed to be in canonical ordering
+/// @returns The ket index of the provided group
+///
+/// @tparam Group The type of the container used to represent the index group.
+/// May be a std::pair<Index, Index> or a container that retains
+/// the order of its elements based on how they have been inserted
+/// (vector-like).
+template <IndexGroup Group>
+decltype(auto) get_ket_idx(Group&& group) {
+  if constexpr (IndexTuple<Group>) {
+    return std::get<1>(group);
+  } else {
+    using std::begin;
+    using std::end;
+    auto it = begin(group);
+    SEQUANT_ASSERT(it != end(group));
     std::advance(it, 1);
+    SEQUANT_ASSERT(it != end(group));
     return *it;
   }
 }
 
+/// Produces a view of plain Index objects from a range of SlottedIndex objects
+///
+/// @note This function will bring indices in the group into canonical order
+/// (based on slot type) or assert that they already are in canonical order,
+/// if group is const.
+decltype(auto) as_index_group_view(SlottedIndexGroup auto&& group) {
+  static_assert(static_cast<int>(SlotType::Bra) == 0);
+  static_assert(static_cast<int>(SlotType::Ket) == 1);
+
+  // We have to ensure a unique order of indices if we're getting rid of the
+  // SlotType tag
+  if constexpr (!meta::is_immutable_v<decltype(group)>) {
+    std::ranges::sort(group, std::less{}, &SlottedIndex::slot_type);
+  } else {
+    SEQUANT_ASSERT(
+        std::ranges::is_sorted(group, std::less{}, &SlottedIndex::slot_type));
+  }
+
+  return group | std::ranges::views::transform(
+                     [](auto&& idx) -> decltype(auto) { return idx.index(); });
+}
+
+/// Produces a view of a view of plain Index objects from a range of a range of
+/// SlottedIndex objects
+///
+/// @note This function transforms all elements via as_index_group, which may
+/// change order of indices or assert on the order being as expected.
+///
+/// @sa as_index_group
+decltype(auto) as_view_of_index_groups(
+    SlottedIndexGroupContainer auto&& groups) {
+  return groups |
+         std::ranges::views::transform([](auto&& group) -> decltype(auto) {
+           return as_index_group_view(group);
+         });
+}
+
 }  // namespace sequant
 
 #endif  // SEQUANT_CORE_UTILITY_INDICES_HPP
diff --git a/SeQuant/core/utility/overloads.hpp b/SeQuant/core/utility/overloads.hpp
new file mode 100644
index 0000000000..60e6d58957
--- /dev/null
+++ b/SeQuant/core/utility/overloads.hpp
@@ -0,0 +1,22 @@
+#ifndef SEQUANT_CORE_UTILITY_OVERLOADS_HPP
+#define SEQUANT_CORE_UTILITY_OVERLOADS_HPP
+
+#include <initializer_list>
+
+namespace sequant {
+
+/// A class that can't be instantiated due to the private constructor.
+/// Useful for enforcing initializer_lists can only ever be empty
+class Nonconstructible {
+ private:
+  Nonconstructible() = default;
+};
+
+/// Tag for providing unambiguous overloads for catching function(bla, {})
+/// where the initializer list could otherwise be matched to multiple other
+/// overloads.
+using EmptyInitializerList = std::initializer_list<Nonconstructible>;
+
+}  // namespace sequant
+
+#endif  // SEQUANT_CORE_UTILITY_OVERLOADS_HPP
diff --git a/SeQuant/domain/mbpt/biorthogonalization.cpp b/SeQuant/domain/mbpt/biorthogonalization.cpp
index 15ff7ee1b9..480bfea2ba 100644
--- a/SeQuant/domain/mbpt/biorthogonalization.cpp
+++ b/SeQuant/domain/mbpt/biorthogonalization.cpp
@@ -1,4 +1,5 @@
 #include <SeQuant/domain/mbpt/biorthogonalization.hpp>
+#include <SeQuant/domain/mbpt/detail/concepts.hpp>
 #include <SeQuant/domain/mbpt/spin.hpp>
 
 #include <SeQuant/core/container.hpp>
@@ -531,11 +532,10 @@ void biorthogonal_transform(container::svector<ResultExpr>& result_exprs,
   }
 }
 
-ExprPtr biorthogonal_transform(
-    const sequant::ExprPtr& expr,
-    const container::svector<container::svector<sequant::Index>>&
-        ext_index_groups,
-    const double threshold) {
+template <detail::index_group_range IdxGroups>
+ExprPtr biorthogonal_transform_impl(const sequant::ExprPtr& expr,
+                                    IdxGroups&& ext_index_groups,
+                                    const double threshold) {
   ResultExpr res(
       bra(ext_index_groups | ranges::views::transform([](const auto& pair) {
             return get_ket_idx(pair);
@@ -553,9 +553,25 @@ ExprPtr biorthogonal_transform(
   return res.expression();
 }
 
-ExprPtr WK_biorthogonalization_filter(
-    ExprPtr expr,
-    const container::svector<container::svector<Index>>& ext_idxs) {
+ExprPtr biorthogonal_transform(
+    const sequant::ExprPtr& expr,
+    const container::svector<container::svector<sequant::SlottedIndex>>&
+        ext_index_groups,
+    const double threshold) {
+  return biorthogonal_transform_impl(
+      expr, as_view_of_index_groups(ext_index_groups), threshold);
+}
+
+ExprPtr biorthogonal_transform(
+    const sequant::ExprPtr& expr,
+    const container::svector<container::svector<sequant::Index>>&
+        ext_index_groups,
+    const double threshold) {
+  return biorthogonal_transform_impl(expr, ext_index_groups, threshold);
+}
+
+template <detail::index_group_range IdxGroups>
+ExprPtr WK_biorthogonalization_filter_impl(ExprPtr expr, IdxGroups&& ext_idxs) {
   if (!expr->is<Sum>()) return expr;
   if (ext_idxs.size() <= 2) return expr;  // always skip R1 and R2
 
@@ -603,10 +619,21 @@ ExprPtr WK_biorthogonalization_filter(
   return result;
 }
 
-ExprPtr biorthogonal_transform_pre_nnsproject(
-    ExprPtr& expr,
-    const container::svector<container::svector<Index>>& ext_idxs,
-    bool factor_out_nns_projector) {
+ExprPtr WK_biorthogonalization_filter(
+    ExprPtr expr,
+    const container::svector<container::svector<SlottedIndex>>& ext_idxs) {
+  return WK_biorthogonalization_filter_impl(expr,
+                                            as_view_of_index_groups(ext_idxs));
+}
+ExprPtr WK_biorthogonalization_filter(
+    ExprPtr expr,
+    const container::svector<container::svector<Index>>& ext_idxs) {
+  return WK_biorthogonalization_filter_impl(expr, ext_idxs);
+}
+
+template <detail::index_group_range IdxGroups>
+ExprPtr biorthogonal_transform_pre_nnsproject_impl(
+    ExprPtr& expr, IdxGroups&& ext_idxs, bool factor_out_nns_projector) {
   using ranges::views::transform;
 
   // Remove leading S operator if present
@@ -616,7 +643,8 @@ ExprPtr biorthogonal_transform_pre_nnsproject(
           remove_tensor(term.as_shared_ptr<Product>(), reserved::symm_label());
   }
 
-  auto bt = biorthogonal_transform(expr, ext_idxs);
+  auto bt = biorthogonal_transform_impl(
+      expr, ext_idxs, default_biorthogonalizer_pseudoinverse_threshold);
 
   auto bixs = ext_idxs | transform([](auto&& vec) { return get_bra_idx(vec); });
   auto kixs = ext_idxs | transform([](auto&& vec) { return get_ket_idx(vec); });
@@ -631,7 +659,7 @@ ExprPtr biorthogonal_transform_pre_nnsproject(
 
     bt = S_maps(bt);
     canonicalize(bt);
-    bt = WK_biorthogonalization_filter(bt, ext_idxs);
+    bt = WK_biorthogonalization_filter_impl(bt, ext_idxs);
   }
 
   bt = S_tensor * bt;
@@ -640,6 +668,22 @@ ExprPtr biorthogonal_transform_pre_nnsproject(
   return bt;
 }
 
+ExprPtr biorthogonal_transform_pre_nnsproject(
+    ExprPtr& expr,
+    const container::svector<container::svector<SlottedIndex>>& ext_idxs,
+    bool factor_out_nns_projector) {
+  return biorthogonal_transform_pre_nnsproject_impl(
+      expr, as_view_of_index_groups(ext_idxs), factor_out_nns_projector);
+}
+
+ExprPtr biorthogonal_transform_pre_nnsproject(
+    ExprPtr& expr,
+    const container::svector<container::svector<Index>>& ext_idxs,
+    bool factor_out_nns_projector) {
+  return biorthogonal_transform_pre_nnsproject_impl(expr, ext_idxs,
+                                                    factor_out_nns_projector);
+}
+
 namespace detail {
 
 std::vector<double> compute_nns_p_coeffs(std::size_t n_particles,
diff --git a/SeQuant/domain/mbpt/biorthogonalization.hpp b/SeQuant/domain/mbpt/biorthogonalization.hpp
index 5593b2d5f1..afb7912475 100644
--- a/SeQuant/domain/mbpt/biorthogonalization.hpp
+++ b/SeQuant/domain/mbpt/biorthogonalization.hpp
@@ -4,6 +4,8 @@
 #include <SeQuant/core/container.hpp>
 #include <SeQuant/core/expr.hpp>
 #include <SeQuant/core/index.hpp>
+#include <SeQuant/core/slotted_index.hpp>
+#include <SeQuant/core/utility/indices.hpp>
 #include <SeQuant/core/utility/memoize.hpp>
 
 #if defined(SEQUANT_HAS_TILEDARRAY)
@@ -43,6 +45,12 @@ void biorthogonal_transform(
 /// @note uses hardcoded coefficients for ranks 1-5,
 ///  for higher ranks computes coefficients (if Eigen3 is available, else throws
 ///  an exception)
+[[nodiscard]] ExprPtr biorthogonal_transform(
+    const ExprPtr& expr,
+    const container::svector<container::svector<sequant::SlottedIndex>>&
+        ext_index_groups = {},
+    double pseudoinverse_threshold =
+        default_biorthogonalizer_pseudoinverse_threshold);
 [[nodiscard]] ExprPtr biorthogonal_transform(
     const ExprPtr& expr,
     const container::svector<container::svector<sequant::Index>>&
@@ -66,6 +74,9 @@ void biorthogonal_transform(
 /// @param ext_idxs A vector of external index groups. The function will not
 /// apply the filtering logic if `ext_idxs.size()` is 2 or less.
 /// @return A new `ExprPtr` representing the filtered and compacted expression.
+ExprPtr WK_biorthogonalization_filter(
+    ExprPtr expr,
+    const container::svector<container::svector<SlottedIndex>>& ext_idxs);
 ExprPtr WK_biorthogonalization_filter(
     ExprPtr expr,
     const container::svector<container::svector<Index>>& ext_idxs);
@@ -83,6 +94,10 @@ ExprPtr WK_biorthogonalization_filter(
 /// projector. If false, NNS projector is not factored out.
 /// @return Expression pointer to the biorthogonalized result with leading S
 /// operator.
+ExprPtr biorthogonal_transform_pre_nnsproject(
+    ExprPtr& expr,
+    const container::svector<container::svector<SlottedIndex>>& ext_idxs,
+    bool factor_out_nns_projector = true);
 ExprPtr biorthogonal_transform_pre_nnsproject(
     ExprPtr& expr,
     const container::svector<container::svector<Index>>& ext_idxs,
diff --git a/SeQuant/domain/mbpt/detail/concepts.hpp b/SeQuant/domain/mbpt/detail/concepts.hpp
new file mode 100644
index 0000000000..7853650512
--- /dev/null
+++ b/SeQuant/domain/mbpt/detail/concepts.hpp
@@ -0,0 +1,20 @@
+#ifndef SEQUANT_DOMAIN_MBPT_DETAIL_CONCEPTS_HPP
+#define SEQUANT_DOMAIN_MBPT_DETAIL_CONCEPTS_HPP
+
+#include <SeQuant/core/index.hpp>
+
+#include <ranges>
+
+namespace sequant::mbpt::detail {
+
+/// Concept for a range of a range of Index objects
+template <typename T>
+concept index_group_range =
+    std::ranges::range<T> &&
+    std::ranges::range<std::ranges::range_value_t<T>> &&
+    std::convertible_to<
+        std::ranges::range_value_t<std::ranges::range_value_t<T>>, Index>;
+
+}  // namespace sequant::mbpt::detail
+
+#endif  // SEQUANT_DOMAIN_MBPT_DETAIL_CONCEPTS_HPP
diff --git a/SeQuant/domain/mbpt/op.hpp b/SeQuant/domain/mbpt/op.hpp
index ed79d8eef8..dbea3ff144 100644
--- a/SeQuant/domain/mbpt/op.hpp
+++ b/SeQuant/domain/mbpt/op.hpp
@@ -109,7 +109,7 @@ struct OpParams {
     if (!batch_ordinals.empty()) {
       SEQUANT_ASSERT(ranges::is_sorted(batch_ordinals) &&
                      "OpParams: batch_ordinals must be sorted");
-      auto duplicate = ranges::adjacent_find(batch_ordinals);
+      [[maybe_unused]] auto duplicate = ranges::adjacent_find(batch_ordinals);
       SEQUANT_ASSERT(duplicate == batch_ordinals.end() &&
                      "OpParams: batch_ordinals must contain unique values");
     }
diff --git a/SeQuant/domain/mbpt/spin.cpp b/SeQuant/domain/mbpt/spin.cpp
index c1c6921ee5..edce56d919 100644
--- a/SeQuant/domain/mbpt/spin.cpp
+++ b/SeQuant/domain/mbpt/spin.cpp
@@ -1,18 +1,22 @@
 #include <SeQuant/domain/mbpt/biorthogonalization.hpp>
+#include <SeQuant/domain/mbpt/detail/concepts.hpp>
 #include <SeQuant/domain/mbpt/spin.hpp>
 
 #include <SeQuant/core/algorithm.hpp>
 #include <SeQuant/core/attr.hpp>
 #include <SeQuant/core/expr.hpp>
+#include <SeQuant/core/index.hpp>
 #include <SeQuant/core/math.hpp>
 #include <SeQuant/core/optimize.hpp>
 #include <SeQuant/core/rational.hpp>
 #include <SeQuant/core/reserved.hpp>
+#include <SeQuant/core/slotted_index.hpp>
 #include <SeQuant/core/space.hpp>
 #include <SeQuant/core/tensor_canonicalizer.hpp>
 #include <SeQuant/core/tensor_network.hpp>
 #include <SeQuant/core/utility/indices.hpp>
 #include <SeQuant/core/utility/macros.hpp>
+#include <SeQuant/core/utility/overloads.hpp>
 #include <SeQuant/core/utility/permutation.hpp>
 #include <SeQuant/core/utility/swap.hpp>
 
@@ -37,7 +41,6 @@
 #include <functional>
 #include <iterator>
 #include <memory>
-#include <new>
 #include <numeric>
 #include <stdexcept>
 #include <string_view>
@@ -107,10 +110,10 @@ template <typename Container, typename TraceFunction, typename... Args>
 
   if (!searchForNonEquivalentResults) {
     ResultExpr traced = expr.clone();
-    traced.expression() =
-        tracer(traced.expression(),
-               traced.index_particle_grouping<container::svector<Index>>(),
-               std::forward<Args>(args)...);
+    traced.expression() = tracer(
+        traced.expression(),
+        traced.index_particle_grouping<container::svector<SlottedIndex>>(),
+        std::forward<Args>(args)...);
 
     traced.set_symmetry(Symmetry::Nonsymm);
 
@@ -198,10 +201,10 @@ template <typename Container, typename TraceFunction, typename... Args>
       }
     }();
 
-    result.expression() =
-        tracer(result.expression(),
-               result.index_particle_grouping<container::svector<Index>>(),
-               std::forward<Args>(args)...);
+    result.expression() = tracer(
+        result.expression(),
+        result.index_particle_grouping<container::svector<SlottedIndex>>(),
+        std::forward<Args>(args)...);
 
     result.set_symmetry(Symmetry::Nonsymm);
 
@@ -926,10 +929,10 @@ ExprPtr S_maps(const ExprPtr& expr) {
   return result;
 }
 
-ExprPtr closed_shell_spintrace(
-    const ExprPtr& expression,
-    const container::svector<container::svector<Index>>& ext_index_groups,
-    bool full_expansion) {
+template <detail::index_group_range IdxGroups>
+ExprPtr closed_shell_spintrace_impl(const ExprPtr& expression,
+                                    IdxGroups&& ext_index_groups,
+                                    bool full_expansion) {
   // Symmetrize and expression
   // Partially expand the antisymmetrizer and write it in terms of S operator.
   // See symmetrize_expr(expr) function for implementation details. We want an
@@ -951,7 +954,7 @@ ExprPtr closed_shell_spintrace(
     rapid_simplify(temp);
     return temp;
   };
-  auto expr = partially_or_fully_expand(expression);
+  ExprPtr expr = partially_or_fully_expand(expression);
 
   // Index tags are cleaned prior to calling the fast canonicalizer
   detail::reset_idx_tags(expr);  // This call is REQUIRED
@@ -1019,11 +1022,11 @@ ExprPtr closed_shell_spintrace(
     };
     auto product_bras = get_bra_indices(temp_product);
 
-    auto substitute_ext_idx = [&product_bras, &product_kets](
-                                  const container::svector<Index>& idx_pair) {
+    auto substitute_ext_idx = [&product_bras,
+                               &product_kets](const auto& idx_pair) {
       SEQUANT_ASSERT(idx_pair.size() == 2);
-      const auto& what = idx_pair[0];
-      const auto& with = idx_pair[1];
+      const auto& what = get_bra_idx(idx_pair);
+      const auto& with = get_ket_idx(idx_pair);
       std::replace(product_bras.begin(), product_bras.end(), what, with);
       std::replace(product_kets.begin(), product_kets.end(), what, with);
     };
@@ -1065,11 +1068,35 @@ ExprPtr closed_shell_spintrace(
   }
 }
 
+ExprPtr closed_shell_spintrace(
+    const ExprPtr& expression,
+    const container::svector<container::svector<SlottedIndex>>&
+        ext_index_groups,
+    bool full_expansion) {
+  return closed_shell_spintrace_impl(
+      expression, as_view_of_index_groups(ext_index_groups), full_expansion);
+}
+
+ExprPtr closed_shell_spintrace(const ExprPtr& expression, EmptyInitializerList,
+                               bool full_expansion) {
+  return closed_shell_spintrace_impl(
+      expression, container::svector<container::svector<Index>>{},
+      full_expansion);
+}
+
+ExprPtr closed_shell_spintrace(
+    const ExprPtr& expression,
+    const container::svector<container::svector<Index>>& ext_index_groups,
+    bool full_expansion) {
+  return closed_shell_spintrace_impl(expression, ext_index_groups,
+                                     full_expansion);
+}
+
 container::svector<ResultExpr> closed_shell_spintrace(const ResultExpr& expr,
                                                       bool full_expansion) {
-  using TraceFunction =
-      ExprPtr (*)(const ExprPtr&,
-                  const container::svector<container::svector<Index>>&, bool);
+  using TraceFunction = ExprPtr (*)(
+      const ExprPtr&,
+      const container::svector<container::svector<SlottedIndex>>&, bool);
 
   return detail::wrap_trace<container::svector<ResultExpr>>(
       expr, static_cast<TraceFunction>(&closed_shell_spintrace),
@@ -1343,10 +1370,10 @@ std::vector<ExprPtr> open_shell_P_op_vector(const Tensor& A) {
   return result_vector;
 }
 
-std::vector<ExprPtr> open_shell_spintrace(
-    const ExprPtr& expr,
-    const container::svector<container::svector<Index>>& ext_index_groups,
-    std::optional<int> target_spin_case) {
+template <detail::index_group_range IdxGroups>
+std::vector<ExprPtr> open_shell_spintrace_impl(
+    const ExprPtr& expr, IdxGroups&& ext_index_groups,
+    const std::optional<int>& target_spin_case) {
   if (expr->is<Constant>() || expr->is<Variable>()) {
     return std::vector<ExprPtr>{expr};
   }
@@ -1356,10 +1383,11 @@ std::vector<ExprPtr> open_shell_spintrace(
       get_used_indices<decltype(grand_idxlist)>(expr);
 
   container::set<Index> ext_idxlist;
-  for (auto&& idxgrp : ext_index_groups) {
-    for (auto&& idx : idxgrp) {
+  for (const auto& idxgrp : ext_index_groups) {
+    for (const Index& current : idxgrp) {
+      Index idx = current;
       idx.reset_tag();
-      ext_idxlist.insert(idx);
+      ext_idxlist.insert(std::move(idx));
     }
   }
 
@@ -1407,26 +1435,25 @@ std::vector<ExprPtr> open_shell_spintrace(
   };
 
   // External index replacement maps
-  auto ext_spin_cases =
-      [&make_spinspecific](const container::svector<IndexGroup>& idx_group) {
-        container::svector<container::map<Index, Index>> all_replacements;
-
-        // container::svector<int> spins(idx_group.size(), 0);
-        for (std::size_t i = 0; i <= idx_group.size(); ++i) {
-          container::svector<int> spins(idx_group.size(), 0);
-          std::fill(spins.end() - i, spins.end(), 1);
-
-          container::map<Index, Index> idx_rep;
-          for (std::size_t j = 0; j != idx_group.size(); ++j) {
-            for (auto& idx : idx_group[j]) {
-              auto spin_idx = make_spinspecific(idx, spins[j]);
-              idx_rep.emplace(idx, spin_idx);
-            }
-          }
-          all_replacements.push_back(idx_rep);
+  auto ext_spin_cases = [&make_spinspecific](const auto& idx_groups) {
+    container::svector<container::map<Index, Index>> all_replacements;
+
+    // container::svector<int> spins(idx_group.size(), 0);
+    for (std::size_t i = 0; i <= idx_groups.size(); ++i) {
+      container::svector<int> spins(idx_groups.size(), 0);
+      std::fill(spins.end() - i, spins.end(), 1);
+
+      container::map<Index, Index> idx_rep;
+      for (std::size_t j = 0; j != idx_groups.size(); ++j) {
+        for (const Index& idx : idx_groups[j]) {
+          auto spin_idx = make_spinspecific(idx, spins[j]);
+          idx_rep.emplace(idx, spin_idx);
         }
-        return all_replacements;
-      };
+      }
+      all_replacements.push_back(idx_rep);
+    }
+    return all_replacements;
+  };
 
   // Internal and external index replacements are independent
   auto i_rep = spin_cases(int_index_groups);
@@ -1485,7 +1512,7 @@ std::vector<ExprPtr> open_shell_spintrace(
     // Loop over internal index replacement maps
     for (auto& i : i_rep) {
       // Add spin labels to internal indices, expand antisymmetric tensors
-      auto spin_expr_i = append_spin(spin_expr, i);
+      ExprPtr spin_expr_i = append_spin(spin_expr, i);
       spin_expr_i = expand_antisymm(spin_expr_i, true);
       expand(spin_expr_i);
       detail::reset_idx_tags(spin_expr_i);
@@ -1529,6 +1556,29 @@ std::vector<ExprPtr> open_shell_spintrace(
   return result;
 }
 
+std::vector<ExprPtr> open_shell_spintrace(
+    const ExprPtr& expr,
+    const container::svector<container::svector<SlottedIndex>>&
+        ext_index_groups,
+    const std::optional<int>& target_spin_case) {
+  return open_shell_spintrace_impl(
+      expr, as_view_of_index_groups(ext_index_groups), target_spin_case);
+}
+
+std::vector<ExprPtr> open_shell_spintrace(
+    const ExprPtr& expr, EmptyInitializerList,
+    const std::optional<int>& target_spin_case) {
+  return open_shell_spintrace_impl(
+      expr, container::svector<container::svector<Index>>{}, target_spin_case);
+}
+
+std::vector<ExprPtr> open_shell_spintrace(
+    const ExprPtr& expr,
+    const container::svector<container::svector<Index>>& ext_index_groups,
+    const std::optional<int>& target_spin_case) {
+  return open_shell_spintrace_impl(expr, ext_index_groups, target_spin_case);
+}
+
 std::vector<ExprPtr> open_shell_CC_spintrace(const ExprPtr& expr) {
   SEQUANT_ASSERT(expr->is<Sum>() || expr->is<Product>());
   Tensor A = expr.is<Sum>() ? expr->at(0)->at(0)->as<Tensor>()
@@ -1579,10 +1629,9 @@ std::vector<ExprPtr> open_shell_CC_spintrace(const ExprPtr& expr) {
   return expr_vec;
 }
 
-ExprPtr spintrace(
-    const ExprPtr& expression,
-    const container::svector<container::svector<Index>>& ext_index_groups,
-    bool spinfree_index_spaces) {
+template <detail::index_group_range IdxGroups>
+ExprPtr spintrace_impl(const ExprPtr& expression, IdxGroups&& ext_index_groups,
+                       bool spinfree_index_spaces) {
   // Escape immediately if expression is a constant
   if (expression->is<Constant>() || expression->is<Variable>()) {
     return expression;
@@ -1599,9 +1648,7 @@ ExprPtr spintrace(
                                [](const auto& list) { return list.size(); });
       return std::accumulate(sizes.begin(), sizes.end(), 0);
     };
-    auto determined_externals =
-        external_indices<container::svector<container::svector<Index>>>(
-            expression);
+    auto determined_externals = external_indices(expression);
 
     SEQUANT_ASSERT(count_indices(ext_index_groups) ==
                    count_indices(determined_externals));
@@ -1668,10 +1715,11 @@ ExprPtr spintrace(
     // List of external indices, i.e. indices that are not summed over Einstein
     // style (indices that are not repeated in an expression)
     container::set<Index> ext_idxlist;
-    for (auto&& idxgrp : ext_index_groups) {
-      for (auto&& idx : idxgrp) {
+    for (const auto& idxgrp : ext_index_groups) {
+      for (const Index& current : idxgrp) {
+        Index idx = current;
         idx.reset_tag();
-        ext_idxlist.insert(idx);
+        ext_idxlist.insert(std::move(idx));
       }
     }
 
@@ -1693,8 +1741,14 @@ ExprPtr spintrace(
     using IndexGroup = container::svector<Index>;
     container::svector<IndexGroup> index_groups;
     for (auto&& i : int_idxlist) index_groups.emplace_back(IndexGroup(1, i));
-    index_groups.insert(index_groups.end(), ext_index_groups.begin(),
-                        ext_index_groups.end());
+    for (const auto& group : ext_index_groups) {
+      IndexGroup target;
+      target.reserve(group.size());
+      for (const Index& current : group) {
+        target.emplace_back(current);
+      }
+      index_groups.emplace_back(std::move(target));
+    }
 
     // EFV: for each spincase (loop over integer from 0 to 2^n-1, n=#of index
     // groups)
@@ -1730,7 +1784,7 @@ ExprPtr spintrace(
         auto st_expr = spintrace_tensor(spin_expr->as<Tensor>());
         result->append(spinfree_index_spaces ? remove_spin(st_expr) : st_expr);
       } else if (spin_expr->is<Product>()) {
-        auto st_expr = spintrace_product(spin_expr.as_shared_ptr<Product>());
+        ExprPtr st_expr = spintrace_product(spin_expr.as_shared_ptr<Product>());
         if (!st_expr->is<Constant>() || st_expr->as<Constant>().value() != 0) {
           result->append(spinfree_index_spaces ? remove_spin(st_expr)
                                                : st_expr);
@@ -1786,13 +1840,35 @@ ExprPtr spintrace(
 
   detail::reset_idx_tags(result);
   return result;
-}  // ExprPtr spintrace
+}
+
+ExprPtr spintrace(const ExprPtr& expression,
+                  const container::svector<container::svector<SlottedIndex>>&
+                      ext_index_groups,
+                  bool spinfree_index_spaces) {
+  return spintrace_impl(expression, as_view_of_index_groups(ext_index_groups),
+                        spinfree_index_spaces);
+}
+
+ExprPtr spintrace(const ExprPtr& expression, EmptyInitializerList,
+                  bool spinfree_index_spaces) {
+  return spintrace_impl(expression,
+                        container::svector<container::svector<Index>>{},
+                        spinfree_index_spaces);
+}
+
+ExprPtr spintrace(
+    const ExprPtr& expression,
+    const container::svector<container::svector<Index>>& ext_index_groups,
+    bool spinfree_index_spaces) {
+  return spintrace_impl(expression, ext_index_groups, spinfree_index_spaces);
+}
 
 container::svector<ResultExpr> spintrace(const ResultExpr& expr,
                                          bool spinfree_index_spaces) {
-  using TraceFunction =
-      ExprPtr (*)(const ExprPtr&,
-                  const container::svector<container::svector<Index>>&, bool);
+  using TraceFunction = ExprPtr (*)(
+      const ExprPtr&,
+      const container::svector<container::svector<SlottedIndex>>&, bool);
 
   return detail::wrap_trace<container::svector<ResultExpr>>(
       expr, static_cast<TraceFunction>(&spintrace), spinfree_index_spaces);
diff --git a/SeQuant/domain/mbpt/spin.hpp b/SeQuant/domain/mbpt/spin.hpp
index 8d10e22b98..fe95b649e3 100644
--- a/SeQuant/domain/mbpt/spin.hpp
+++ b/SeQuant/domain/mbpt/spin.hpp
@@ -12,10 +12,10 @@
 #include <SeQuant/core/container.hpp>
 #include <SeQuant/core/expr.hpp>
 #include <SeQuant/core/index.hpp>
+#include <SeQuant/core/utility/indices.hpp>
 #include <SeQuant/core/utility/macros.hpp>
+#include <SeQuant/core/utility/overloads.hpp>
 
-#include <cstddef>
-#include <initializer_list>
 #include <string>
 #include <vector>
 
@@ -229,7 +229,14 @@ ExprPtr S_maps(const ExprPtr& expr);
 // clang-format on
 ExprPtr closed_shell_spintrace(
     const ExprPtr& expr,
-    const container::svector<container::svector<Index>>& ext_index_groups = {},
+    const container::svector<container::svector<SlottedIndex>>&
+        ext_index_groups = {},
+    bool full_expansion = false);
+ExprPtr closed_shell_spintrace(const ExprPtr& expr, EmptyInitializerList,
+                               bool full_expansion = false);
+ExprPtr closed_shell_spintrace(
+    const ExprPtr& expr,
+    const container::svector<container::svector<Index>>& ext_index_groups,
     bool full_expansion = false);
 
 container::svector<ResultExpr> closed_shell_spintrace(
@@ -334,10 +341,18 @@ std::vector<ExprPtr> open_shell_P_op_vector(const Tensor& A);
 /// @note this performs full expansion of the antisymmetrizer
 /// @sa open_shell_CC_spintrace
 // clang-format on
+std::vector<ExprPtr> open_shell_spintrace(
+    const ExprPtr& expr,
+    const container::svector<container::svector<SlottedIndex>>&
+        ext_index_groups,
+    const std::optional<int>& target_spin_case = std::nullopt);
+std::vector<ExprPtr> open_shell_spintrace(
+    const ExprPtr& expr, EmptyInitializerList,
+    const std::optional<int>& target_spin_case = std::nullopt);
 std::vector<ExprPtr> open_shell_spintrace(
     const ExprPtr& expr,
     const container::svector<container::svector<Index>>& ext_index_groups,
-    std::optional<int> target_spin_case = std::nullopt);
+    const std::optional<int>& target_spin_case = std::nullopt);
 
 // clang-format off
 /// @brief Like open_shell_spintrace but uses minimal expansion of the antisymmetrizer
@@ -370,9 +385,15 @@ std::vector<ExprPtr> open_shell_CC_spintrace(const ExprPtr& expr);
 /// @return an expression with spin integrated/adapted
 /// @warning The result of this function is not simplified since this is a
 /// building block for more specialized spin-tracing functions
+ExprPtr spintrace(const ExprPtr& expr,
+                  const container::svector<container::svector<SlottedIndex>>&
+                      ext_index_groups = {},
+                  bool spinfree_index_spaces = true);
+ExprPtr spintrace(const ExprPtr& expr, EmptyInitializerList,
+                  bool spinfree_index_spaces = true);
 ExprPtr spintrace(
     const ExprPtr& expr,
-    const container::svector<container::svector<Index>>& ext_index_groups = {},
+    const container::svector<container::svector<Index>>& ext_index_groups,
     bool spinfree_index_spaces = true);
 
 container::svector<ResultExpr> spintrace(const ResultExpr& expr,
diff --git a/tests/integration/srcc.cpp b/tests/integration/srcc.cpp
index ee8bdc65f2..a84bdc013e 100644
--- a/tests/integration/srcc.cpp
+++ b/tests/integration/srcc.cpp
@@ -99,8 +99,12 @@ class compute_cceqvec {
         eqvec_sf_ref[R] = closed_shell_spintrace(eqvec_so[R], ext_idxs);
         if (R == 1) {  // closed_shell_spintrace omits 1-body S
           using ranges::views::transform;
-          auto bixs = ext_idxs | transform([](auto&& vec) { return vec[0]; });
-          auto kixs = ext_idxs | transform([](auto&& vec) { return vec[1]; });
+          auto bixs = ext_idxs | transform([](auto&& group) {
+                        return get_bra_idx(group);
+                      });
+          auto kixs = ext_idxs | transform([](auto&& group) {
+                        return get_ket_idx(group);
+                      });
           auto s_tensor =
               ex<Tensor>(Tensor{reserved::symm_label(), bra(kixs), ket(bixs)});
           eqvec_sf_ref[R] = s_tensor * eqvec_sf_ref[R];
diff --git a/tests/unit/test_spin.cpp b/tests/unit/test_spin.cpp
index 5cf8a082c7..7aa7cc8e92 100644
--- a/tests/unit/test_spin.cpp
+++ b/tests/unit/test_spin.cpp
@@ -41,6 +41,8 @@ TEST_CASE("spin", "[spin]") {
   using reserved::antisymm_label;
   using reserved::symm_label;
 
+  using IdxGroupList = container::svector<container::svector<Index>>;
+
   // P.S. ref outputs produced with complete canonicalization
   auto ctx = get_default_context();
   ctx.set(CanonicalizeOptions{.method = CanonicalizationMethod::Complete});
@@ -264,7 +266,8 @@ TEST_CASE("spin", "[spin]") {
       const auto expr = ex<Constant>(rational{1, 4}) *
                         ex<Tensor>(L"g", bra{L"p_1", L"p_2"},
                                    ket{L"p_3", L"p_4"}, Symmetry::Antisymm);
-      auto result = spintrace(expr, {{L"p_1", L"p_3"}, {L"p_2", L"p_4"}});
+      auto result =
+          spintrace(expr, IdxGroupList{{L"p_1", L"p_3"}, {L"p_2", L"p_4"}});
       REQUIRE_THAT(result,
                    EquivalentTo("-1/2 g{p1,p2;p4,p3} + g{p1,p2;p3,p4}"));
     }
@@ -275,7 +278,8 @@ TEST_CASE("spin", "[spin]") {
       // ends up as a single tensor (due to an additional factor of -1 induces
       // by the chosen pairing of external indices)
       ExprPtr expr = parse_expr(L"- g{p1,p2;p3,p4}:A");
-      ExprPtr result = spintrace(expr, {{L"p_1", L"p_4"}, {L"p_2", L"p_3"}});
+      ExprPtr result =
+          spintrace(expr, IdxGroupList{{L"p_1", L"p_4"}, {L"p_2", L"p_3"}});
 
       REQUIRE_THAT(result,
                    EquivalentTo(L"4 g{p1,p2;p4,p3} - 2 g{p1,p2;p3,p4}"));
@@ -284,8 +288,9 @@ TEST_CASE("spin", "[spin]") {
       const auto expr = ex<Constant>(rational{1, 4}) *
                         ex<Tensor>(L"g", bra{L"p_1", L"p_2"},
                                    ket{L"p_3", L"p_4"}, Symmetry::Antisymm);
-      auto result = spintrace(expr, {{L"p_1", L"p_3"}, {L"p_2", L"p_4"}},
-                              /* spinfree = */ false);
+      auto result =
+          spintrace(expr, IdxGroupList{{L"p_1", L"p_3"}, {L"p_2", L"p_4"}},
+                    /* spinfree = */ false);
       REQUIRE_THAT(
           result,
           EquivalentTo(
@@ -609,8 +614,8 @@ SECTION("partial spintracing + S_maps = full spintracing") {
                           ket{L"a_1", L"a_2"}, Symmetry::Antisymm) *
                ex<Tensor>(L"t", bra{L"a_1", L"a_2"}, ket{L"i_1", L"i_2"},
                           Symmetry::Antisymm);
-  auto result =
-      closed_shell_spintrace(input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}});
+  auto result = closed_shell_spintrace(
+      input, IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}});
   REQUIRE_THAT(
       result,
       EquivalentTo("-2 Ŝ{i_1,i_2;a_1,a_2}:N-C-S * t{a_1,a_2;i_2,i_1}:N-C-S "
@@ -621,8 +626,8 @@ SECTION("partial spintracing + S_maps = full spintracing") {
                   "-1 t{a_1,a_2;i_2,i_1}:N-C-S - 1 t{a_2,a_1;i_1,i_2}:N-C-S "
                   "+ 2 t{a_1,a_2;i_1,i_2}:N-C-S + 2 t{a_2,a_1;i_2,i_1}:N-C-S"));
 
-  result =
-      closed_shell_spintrace(input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}}, true);
+  result = closed_shell_spintrace(
+      input, IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}}, true);
   REQUIRE_THAT(
       result, EquivalentTo(
                   "-1 t{a_1,a_2;i_2,i_1}:N-C-S - 1 t{a_2,a_1;i_1,i_2}:N-C-S "
@@ -824,8 +829,8 @@ SECTION("Closed-shell spintrace CCSD") {
     // A * f
     const auto input = ex<Tensor>(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) *
                        ex<Tensor>(L"f", bra{L"a_1"}, ket{L"i_1"});
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     rapid_simplify(result);
     canonicalize(result);
@@ -838,8 +843,8 @@ SECTION("Closed-shell spintrace CCSD") {
                        ex<Tensor>(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) *
                        ex<Tensor>(L"f", bra{L"i_2"}, ket{L"i_1"}) *
                        ex<Tensor>(L"t", bra{L"a_1"}, ket{L"i_2"});
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(result, EquivalentTo("- f{i2;i1} t{a1;i2}"));
   }
@@ -849,8 +854,8 @@ SECTION("Closed-shell spintrace CCSD") {
     const auto input = ex<Tensor>(antisymm_label(), bra{L"i_1"}, ket{L"a_1"}) *
                        ex<Tensor>(L"f", bra{L"a_1"}, ket{L"a_2"}) *
                        ex<Tensor>(L"t", bra{L"a_2"}, ket{L"i_1"});
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(result, EquivalentTo("f{a1;a2} t{a2;i1}"));
   }
@@ -863,8 +868,8 @@ SECTION("Closed-shell spintrace CCSD") {
                                   ket{L"i_1", L"a_2"}, Symmetry::Antisymm) *
                        ex<Tensor>(L"t", bra{L"a_1", L"a_2"},
                                   ket{L"i_2", L"i_3"}, Symmetry::Antisymm);
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(
         result,
@@ -880,8 +885,8 @@ SECTION("Closed-shell spintrace CCSD") {
                                   ket{L"a_2", L"a_3"}, Symmetry::Antisymm) *
                        ex<Tensor>(L"t", bra{L"a_2", L"a_3"},
                                   ket{L"i_1", L"i_2"}, Symmetry::Antisymm);
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
 
     REQUIRE_THAT(result, EquivalentTo("- g{a1,i2;a3,a2} t{a2,a3;i1,i2} + 2 "
@@ -895,8 +900,8 @@ SECTION("Closed-shell spintrace CCSD") {
         ex<Tensor>(L"f", bra{L"i_2"}, ket{L"a_2"}, Symmetry::Antisymm) *
         ex<Tensor>(L"t", bra{L"a_1", L"a_2"}, ket{L"i_1", L"i_2"},
                    Symmetry::Antisymm);
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(
         result,
@@ -910,8 +915,8 @@ SECTION("Closed-shell spintrace CCSD") {
                                   ket{L"a_2", L"a_3"}, Symmetry::Antisymm) *
                        ex<Tensor>(L"t", bra{L"a_2"}, ket{L"i_2"}) *
                        ex<Tensor>(L"t", bra{L"a_3"}, ket{L"i_1"});
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(result, EquivalentTo("2 g{a1,i2;a3,a2} t{a2;i2} t{a3;i1} - "
                                       "g{a1,i2;a3,a2} t{a2;i1} t{a3;i2}"));
@@ -924,8 +929,8 @@ SECTION("Closed-shell spintrace CCSD") {
                                   ket{L"i_1", L"a_2"}, Symmetry::Antisymm) *
                        ex<Tensor>(L"t", bra{L"a_2"}, ket{L"i_2"}) *
                        ex<Tensor>(L"t", bra{L"a_1"}, ket{L"i_3"});
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(result, EquivalentTo("g{i2,i3;i1,a2} t{a1;i3} t{a2;i2} - 2 "
                                       "g{i2,i3;a2,i1} t{a1;i3} t{a2;i2}"));
@@ -938,8 +943,8 @@ SECTION("Closed-shell spintrace CCSD") {
                        ex<Tensor>(L"f", bra{L"i_2"}, ket{L"a_2"}) *
                        ex<Tensor>(L"t", bra{L"a_2"}, ket{L"i_1"}) *
                        ex<Tensor>(L"t", bra{L"a_1"}, ket{L"i_2"});
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(result, EquivalentTo("-f{i2;a2} t{a1;i2} t{a2;i1}"));
   }
@@ -953,8 +958,8 @@ SECTION("Closed-shell spintrace CCSD") {
                        ex<Tensor>(L"t", bra{L"a_1"}, ket{L"i_2"}) *
                        ex<Tensor>(L"t", bra{L"a_2", L"a_3"},
                                   ket{L"i_1", L"i_3"}, Symmetry::Antisymm);
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(result,
                  EquivalentTo("g{i2,i3;a2,a3} t{a1;i3} t{a2,a3;i1,i2} - 2 "
@@ -970,8 +975,8 @@ SECTION("Closed-shell spintrace CCSD") {
                        ex<Tensor>(L"t", bra{L"a_2"}, ket{L"i_1"}) *
                        ex<Tensor>(L"t", bra{L"a_1", L"a_3"},
                                   ket{L"i_2", L"i_3"}, Symmetry::Antisymm);
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(result,
                  EquivalentTo("g{i2,i3;a2,a3} t{a2;i1} t{a1,a3;i3,i2} - 2 "
@@ -987,8 +992,8 @@ SECTION("Closed-shell spintrace CCSD") {
                        ex<Tensor>(L"t", bra{L"a_2"}, ket{L"i_2"}) *
                        ex<Tensor>(L"t", bra{L"a_1", L"a_3"},
                                   ket{L"i_1", L"i_3"}, Symmetry::Antisymm);
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(result,
                  EquivalentTo("-2 g{i2,i3;a2,a3} t{a3;i2} t{a1,a2;i1,i3} + 4 "
@@ -1005,8 +1010,8 @@ SECTION("Closed-shell spintrace CCSD") {
                  ex<Tensor>(L"t", bra{L"a_2"}, ket{L"i_2"}) *
                  ex<Tensor>(L"t", bra{L"a_3"}, ket{L"i_1"}) *
                  ex<Tensor>(L"t", bra{L"a_1"}, ket{L"i_3"});
-    auto result =
-        ex<Constant>(rational{1, 2}) * spintrace(input, {{L"i_1", L"a_1"}});
+    auto result = ex<Constant>(rational{1, 2}) *
+                  spintrace(input, IdxGroupList{{L"i_1", L"a_1"}});
     expand(result);
     REQUIRE_THAT(result,
                  EquivalentTo("-2 g{i2,i3;a2,a3} t{a1;i3} t{a2;i2} t{a3;i1} + "
@@ -1071,7 +1076,8 @@ SECTION("Closed-shell spintrace CCSDT terms") {
     auto result = expand_A_op(input);
     REQUIRE(result->size() == 36);
     result = closed_shell_spintrace(
-        input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
+        input,
+        IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
     REQUIRE(result->size() == 4);
 
     REQUIRE_THAT(
@@ -1101,7 +1107,9 @@ SECTION("Closed-shell spintrace CCSDT terms") {
     REQUIRE(result_1->size() == 36);
     result_1 = expand_antisymm(result_1);
     result_1 = closed_shell_spintrace(
-        input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}}, true);
+        input,
+        IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}},
+        true);
     simplify(result_1);
     REQUIRE(result_1->size() == 18);  // 18 raw terms
     REQUIRE_THAT(
@@ -1131,7 +1139,8 @@ SECTION("Closed-shell spintrace CCSDT terms") {
     // expanding by S_map (this method is used in
     // closed_shell_CC_spintrace_v2)
     auto result_2 = closed_shell_spintrace(
-        input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
+        input,
+        IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
     simplify(result_2);
     result_2 = S_maps(result_2);
     simplify(result_2);
@@ -1204,7 +1213,8 @@ SECTION("Closed-shell spintrace CCSDT terms") {
     REQUIRE(result->size() == 2);
     result = expand_antisymm(result);
     result = closed_shell_spintrace(
-        input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
+        input,
+        IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
     REQUIRE(result->size() == 6);
     simplify(result);
     REQUIRE(result->size() == 6);
@@ -1222,7 +1232,8 @@ SECTION("Closed-shell spintrace CCSDT terms") {
     REQUIRE(result->size() == 36);
     result = expand_antisymm(result);
     result = closed_shell_spintrace(
-        input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
+        input,
+        IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
     REQUIRE(result->size() == 20);
   }
 
@@ -1235,7 +1246,8 @@ SECTION("Closed-shell spintrace CCSDT terms") {
     REQUIRE(result->size() == 2);
     result = expand_antisymm(result);
     result = closed_shell_spintrace(
-        input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
+        input,
+        IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
     REQUIRE(result->size() == 12);
     simplify(result);
     REQUIRE(result->size() == 12);
@@ -1477,8 +1489,8 @@ SECTION("Open-shell spin-tracing") {
     auto input = ex<Constant>(rational{1, 4}) *
                  ex<Tensor>(L"g", bra{L"a_1", L"a_2"}, ket{L"i_1", L"i_2"},
                             Symmetry::Antisymm);
-    auto result =
-        open_shell_spintrace(input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}});
+    auto result = open_shell_spintrace(
+        input, IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}});
     REQUIRE(result.size() == 3);
     REQUIRE_THAT(result[0], EquivalentTo("1/4 g{a↑1,a↑2;i↑1,i↑2}:A"));
     REQUIRE_THAT(result[1], EquivalentTo("1/4 g{a↑1,a↓2;i↑1,i↓2}"));
@@ -1492,8 +1504,8 @@ SECTION("Open-shell spin-tracing") {
                  ex<Tensor>(L"t", bra{L"a_1", L"a_2"}, ket{L"i_2", L"i_3"},
                             Symmetry::Antisymm);
 
-    auto result =
-        open_shell_spintrace(input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}});
+    auto result = open_shell_spintrace(
+        input, IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}});
     REQUIRE(result.size() == 3);
     REQUIRE_THAT(result[0],
                  EquivalentTo("1/2 f{i↑3;i↑1} t{a↑1,a↑2;i↑2,i↑3}:A"));
@@ -1508,8 +1520,8 @@ SECTION("Open-shell spin-tracing") {
                  ex<Tensor>(L"g", bra{L"i_3", L"a_1"}, ket{L"i_1", L"i_2"},
                             Symmetry::Antisymm) *
                  ex<Tensor>(L"t", bra{L"a_2"}, ket{L"i_3"}, Symmetry::Nonsymm);
-    auto result =
-        open_shell_spintrace(input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}});
+    auto result = open_shell_spintrace(
+        input, IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}});
     REQUIRE(result.size() == 3);
     REQUIRE(
         toUtf8(to_latex(result[0])) ==
@@ -1530,7 +1542,8 @@ SECTION("Open-shell spin-tracing") {
                  ex<Tensor>(L"t", bra{L"a_2", L"a_3", L"a_4"},
                             ket{L"i_1", L"i_2", L"i_3"}, Symmetry::Antisymm);
     auto result = open_shell_spintrace(
-        input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
+        input,
+        IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
     REQUIRE(result.size() == 4);
     REQUIRE_THAT(result[0],
                  EquivalentTo("1/12 f{a↑1;a↑4} t{a↑2,a↑3,a↑4;i↑1,i↑2,i↑3}:A"));
@@ -1585,7 +1598,8 @@ SECTION("Open-shell spin-tracing") {
                             ket{L"i_3", L"i_4", L"i_5"}, Symmetry::Antisymm);
 
     auto result = open_shell_spintrace(
-        input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
+        input,
+        IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
     REQUIRE(result[0]->size() == 3);
     auto A3_aaa = Tensor(antisymm_label(), bra{i1A, i2A, i3A},
                          ket{a1A, a2A, a3A}, Symmetry::Antisymm);
@@ -1633,7 +1647,8 @@ SECTION("Open-shell spin-tracing") {
     input = expand_P_op(input);
     input->visit(reset_idx_tags);
     auto result = open_shell_spintrace(
-        input, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
+        input,
+        IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
 
     auto result_aab = ex<Tensor>(Tensor(antisymm_label(), bra{i1A, i2A},
                                         ket{a1A, a2A}, Symmetry::Antisymm)) *
@@ -1656,7 +1671,8 @@ SECTION("Open-shell spin-tracing") {
                              ket{L"i_3", L"i_4", L"i_5"}, Symmetry::Antisymm);
 
     auto result2 = open_shell_spintrace(
-        input2, {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
+        input2,
+        IdxGroupList{{L"i_1", L"a_1"}, {L"i_2", L"a_2"}, {L"i_3", L"a_3"}});
     REQUIRE(result2[1]->size() == 24);
   }
 }
diff --git a/tests/unit/test_utilities.cpp b/tests/unit/test_utilities.cpp
index 26cc79a190..2f95d34ea6 100644
--- a/tests/unit/test_utilities.cpp
+++ b/tests/unit/test_utilities.cpp
@@ -8,9 +8,12 @@
 #include <SeQuant/core/index.hpp>
 #include <SeQuant/core/parse.hpp>
 #include <SeQuant/core/reserved.hpp>
+#include <SeQuant/core/slotted_index.hpp>
 #include <SeQuant/core/utility/conversion.hpp>
+#include <SeQuant/core/utility/exception.hpp>
 #include <SeQuant/core/utility/expr.hpp>
 #include <SeQuant/core/utility/indices.hpp>
+#include <SeQuant/core/utility/macros.hpp>
 #include <SeQuant/core/utility/singleton.hpp>
 #include <SeQuant/core/utility/strong.hpp>
 #include <SeQuant/core/utility/tensor.hpp>
@@ -536,41 +539,62 @@ TEST_CASE("utilities", "[utilities]") {
   }
 
   SECTION("external_indices") {
-    for (const auto& [input, expected] :
-         std::vector<std::pair<std::wstring, std::vector<std::vector<Index>>>>{
+    for (const auto& [input, expected] : std::vector<
+             std::pair<std::wstring, std::vector<std::vector<SlottedIndex>>>>{
              {L"1/2", {}},
              {L"t", {}},
              {L"t{}", {}},
-             {L"t{i1;a1}", {{L"i_1", L"a_1"}}},
-             {L"t{i1,i2;a1,a2}", {{L"i_1", L"a_1"}, {L"i_2", L"a_2"}}},
-             {L"t{i1,i2;a2,a1}", {{L"i_1", L"a_2"}, {L"i_2", L"a_1"}}},
-             {L"t{;;i1,i2}", {{L"i_2"}, {L"i_1"}}},
-             {L"t{a1;a2;i1,i2}", {{L"a_1", L"a_2"}, {L"i_1"}, {L"i_2"}}},
+             {L"t{i1;a1}",
+              {{{L"i_1", SlotType::Bra}, {L"a_1", SlotType::Ket}}}},
+             {L"t{i1,i2;a1,a2}",
+              {{{L"i_1", SlotType::Bra}, {L"a_1", SlotType::Ket}},
+               {{L"i_2", SlotType::Bra}, {L"a_2", SlotType::Ket}}}},
+             {L"t{i1,i2;a2,a1}",
+              {{{L"i_1", SlotType::Bra}, {L"a_2", SlotType::Ket}},
+               {{L"i_2", SlotType::Bra}, {L"a_1", SlotType::Ket}}}},
+             {L"t{;;i1,i2}",
+              {{{L"i_2", SlotType::Aux}}, {{L"i_1", SlotType::Aux}}}},
+             {L"t{a1;a2;i1,i2}",
+              {{{L"a_1", SlotType::Bra}, {L"a_2", SlotType::Ket}},
+               {{L"i_1", SlotType::Aux}},
+               {{L"i_2", SlotType::Aux}}}},
              {L"t{a1,a2;a3,a4;i1,i2}",
-              {{L"a_1", L"a_3"}, {L"a_2", L"a_4"}, {L"i_1"}, {L"i_2"}}},
+              {{{L"a_1", SlotType::Bra}, {L"a_3", SlotType::Ket}},
+               {{L"a_2", SlotType::Bra}, {L"a_4", SlotType::Ket}},
+               {{L"i_1", SlotType::Aux}},
+               {{L"i_2", SlotType::Aux}}}},
              {L"g{a3;a2;i2} t{a1,a2;a3,a4;i1,i2}",
-              {{L"a_1", L"a_4"}, {L"i_1"}}},
+              {{{L"a_1", SlotType::Bra}, {L"a_4", SlotType::Ket}},
+               {{L"i_1", SlotType::Aux}}}},
              {L"t{a1,a2;i1,i2} + t{a1;i1} t{a2;i2}",
-              {{L"a_1", L"i_1"}, {L"a_2", L"i_2"}}},
+              {{{L"a_1", SlotType::Bra}, {L"i_1", SlotType::Ket}},
+               {{L"a_2", SlotType::Bra}, {L"i_2", SlotType::Ket}}}},
              // When present, external indices are deduced from the symmetrizer
              // (though bra/ket will be swapped)
-             {L"Â{a1;i1} t{i1;a1}", {{L"i_1", L"a_1"}}},
+             {L"Â{a1;i1} t{i1;a1}",
+              {{{L"i_1", SlotType::Bra}, {L"a_1", SlotType::Ket}}}},
              {L"Â{a1,a2;i1,i2} t{i1;a1} t{i2;a2}",
-              {{L"i_1", L"a_1"}, {L"i_2", "a_2"}}},
-             {L"Ŝ{a1;i1} t{i1;a1}", {{L"i_1", L"a_1"}}},
+              {{{L"i_1", SlotType::Bra}, {L"a_1", SlotType::Ket}},
+               {{L"i_2", SlotType::Bra}, {L"a_2", SlotType::Ket}}}},
+             {L"Ŝ{a1;i1} t{i1;a1}",
+              {{{L"i_1", SlotType::Bra}, {L"a_1", SlotType::Ket}}}},
              {L"Ŝ{a1,a2;i1,i2} t{i1;a1} t{i2;a2}",
-              {{L"i_1", L"a_1"}, {L"i_2", "a_2"}}},
+              {{{L"i_1", SlotType::Bra}, {L"a_1", SlotType::Ket}},
+               {{L"i_2", SlotType::Bra}, {L"a_2", SlotType::Ket}}}},
              // We still want the "swap" behavior when called on an isolated
              // symmetrizer Note: this is the inverse behavior to a regular
              // tensor
-             {L"Ŝ{a1,a2;i1,i2}", {{L"i_1", L"a_1"}, {L"i_2", "a_2"}}},
-             {L"Â{a1;i1}", {{L"i_1", L"a_1"}}},
+             {L"Ŝ{a1,a2;i1,i2}",
+              {{{L"i_1", SlotType::Bra}, {L"a_1", SlotType::Ket}},
+               {{L"i_2", SlotType::Bra}, {L"a_2", SlotType::Ket}}}},
+             {L"Â{a1;i1}",
+              {{{L"i_1", SlotType::Bra}, {L"a_1", SlotType::Ket}}}},
          }) {
       CAPTURE(toUtf8(input));
 
       ExprPtr expr = parse_expr(input);
-      auto actual =
-          external_indices<std::remove_cvref_t<decltype(expected)>>(expr);
+      SlottedIndexGroupContainer auto actual =
+          external_indices<std::vector, std::vector>(expr);
 
       REQUIRE_THAT(actual, ::Catch::Matchers::UnorderedRangeEquals(expected));
     }
@@ -632,4 +656,113 @@ TEST_CASE("utilities", "[utilities]") {
       }
     }
   }
+
+  SECTION("as_index_group_view") {
+    std::vector<SlottedIndex> indices;
+
+    SECTION("as-is") {
+      SECTION("mutable") {
+        indices = {{L"a_1", SlotType::Bra}, {L"i_1", SlotType::Ket}};
+
+        auto group_view = as_index_group_view(indices);
+
+        REQUIRE(get_bra_idx(group_view) == Index(L"a_1"));
+        REQUIRE(get_ket_idx(group_view) == Index(L"i_1"));
+      }
+      SECTION("const") {
+        indices = {{L"a_1", SlotType::Bra}, {L"i_1", SlotType::Ket}};
+
+        auto group_view = as_index_group_view(std::as_const(indices));
+
+        REQUIRE(get_bra_idx(group_view) == Index(L"a_1"));
+        REQUIRE(get_ket_idx(group_view) == Index(L"i_1"));
+      }
+    }
+    SECTION("reordering") {
+      SECTION("mutable") {
+        indices = {{L"a_1", SlotType::Ket}, {L"i_1", SlotType::Bra}};
+
+        auto group_view = as_index_group_view(indices);
+
+        REQUIRE(get_bra_idx(group_view) == Index(L"i_1"));
+        REQUIRE(get_ket_idx(group_view) == Index(L"a_1"));
+      }
+#if SEQUANT_ASSERT_BEHAVIOR == SEQUANT_ASSERT_THROW
+      SECTION("const") {
+        indices = {{L"a_1", SlotType::Ket}, {L"i_1", SlotType::Bra}};
+
+        // Since the indices are not in the expected order but the range is
+        // passed as const, the assertions will trigger due to uncorrectable
+        // ordering.
+        REQUIRE_THROWS_AS(as_index_group_view(std::as_const(indices)),
+                          Exception);
+      }
+#endif
+    }
+  }
+
+  SECTION("as_view_of_index_groups") {
+    std::vector<container::svector<SlottedIndex>> groups;
+
+    SECTION("as-is") {
+      SECTION("mutable") {
+        groups = {
+            {{L"a_1", SlotType::Bra}, {L"i_1", SlotType::Ket}},
+            {{L"a_2", SlotType::Bra}, {L"i_2", SlotType::Ket}},
+        };
+
+        auto groups_view = as_view_of_index_groups(groups);
+
+        REQUIRE(get_bra_idx(groups_view[0]) == Index(L"a_1"));
+        REQUIRE(get_ket_idx(groups_view[0]) == Index(L"i_1"));
+        REQUIRE(get_bra_idx(groups_view[1]) == Index(L"a_2"));
+        REQUIRE(get_ket_idx(groups_view[1]) == Index(L"i_2"));
+      }
+      SECTION("const") {
+        groups = {
+            {{L"a_1", SlotType::Bra}, {L"i_1", SlotType::Ket}},
+            {{L"a_2", SlotType::Bra}, {L"i_2", SlotType::Ket}},
+        };
+
+        auto groups_view = as_view_of_index_groups(std::as_const(groups));
+
+        REQUIRE(get_bra_idx(groups_view[0]) == Index(L"a_1"));
+        REQUIRE(get_ket_idx(groups_view[0]) == Index(L"i_1"));
+        REQUIRE(get_bra_idx(groups_view[1]) == Index(L"a_2"));
+        REQUIRE(get_ket_idx(groups_view[1]) == Index(L"i_2"));
+      }
+    }
+    SECTION("reordering") {
+      SECTION("mutable") {
+        groups = {
+            {{L"a_1", SlotType::Bra}, {L"i_1", SlotType::Ket}},
+            {{L"a_2", SlotType::Ket}, {L"i_2", SlotType::Bra}},
+        };
+
+        auto groups_view = as_view_of_index_groups(groups);
+
+        REQUIRE(get_bra_idx(groups_view[0]) == Index(L"a_1"));
+        REQUIRE(get_ket_idx(groups_view[0]) == Index(L"i_1"));
+        REQUIRE(get_bra_idx(groups_view[1]) == Index(L"i_2"));
+        REQUIRE(get_ket_idx(groups_view[1]) == Index(L"a_2"));
+      }
+      SECTION("const") {
+        groups = {
+            {{L"a_1", SlotType::Bra}, {L"i_1", SlotType::Ket}},
+            {{L"a_2", SlotType::Ket}, {L"i_2", SlotType::Bra}},
+        };
+
+        auto groups_view = as_view_of_index_groups(std::as_const(groups));
+
+        // Note: due to the lazy evaluation of ranges, we get the assert error
+        // only once we actually attempt to access the second group
+        REQUIRE(get_bra_idx(groups_view[0]) == Index(L"a_1"));
+        REQUIRE(get_ket_idx(groups_view[0]) == Index(L"i_1"));
+
+#if SEQUANT_ASSERT_BEHAVIOR == SEQUANT_ASSERT_THROW
+        REQUIRE_THROWS_AS(get_bra_idx(groups_view[1]), Exception);
+#endif
+      }
+    }
+  }
 }