diff --git a/SeQuant/core/utility/indices.hpp b/SeQuant/core/utility/indices.hpp
index 3d090c5d09..31c10a2eae 100644
--- a/SeQuant/core/utility/indices.hpp
+++ b/SeQuant/core/utility/indices.hpp
@@ -103,6 +103,12 @@ IndexGroups<Container> get_unique_indices(const Variable&) {
   return {};
 }
 
+template <typename Container = std::vector<Index>>
+IndexGroups<Container> get_unique_indices([[maybe_unused]] const Power& power) {
+  SEQUANT_ASSERT(power.base()->is<Constant>() || power.base()->is<Variable>());
+  return {};
+}
+
 /// @returns Lists of non-contracted indices arising when contracting the two
 /// given tensors in the order bra, ket, auxiliary
 template <typename Container = container::svector<Index>>
@@ -249,6 +255,8 @@ IndexGroups<Container> get_unique_indices(const Expr& expr) {
     return get_unique_indices<Container>(expr.as<Sum>());
   } else if (expr.is<Product>()) {
     return get_unique_indices<Container>(expr.as<Product>());
+  } else if (expr.is<Power>()) {
+    return get_unique_indices<Container>(expr.as<Power>());
   } else {
     throw Exception(
         "Encountered unsupported expression type in get_unique_indices");
diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp
index 382ece977f..8c61271eba 100644
--- a/SeQuant/domain/mbpt/models/cc.cpp
+++ b/SeQuant/domain/mbpt/models/cc.cpp
@@ -322,6 +322,11 @@ std::vector<ExprPtr> CC::λʼ(size_t rank, size_t order,
   return result;
 }
 
+namespace {
+// EOM eigenvector operators R and L use SquareRoot normalization s
+constexpr Normalization eom_norm = Normalization::SquareRoot;
+}  // namespace
+
 std::vector<ExprPtr> CC::eom_r(nₚ np, nₕ nh) const {
   SEQUANT_ASSERT((np > 0 || nh > 0) && "Unsupported excitation order");
   if (np != nh)
@@ -341,9 +346,9 @@ std::vector<ExprPtr> CC::eom_r(nₚ np, nₕ nh) const {
   // for unitary ansatz, we need to compute the commutator [hbar, R], otherwise
   // just hbar * R is sufficient because ref_av uses connectivity
   if (this->unitary()) {
-    hbar_R = commutator(hbar, R(np, nh));
+    hbar_R = commutator(hbar, R(np, nh, eom_norm));
   } else {
-    hbar_R = hbar * R(np, nh);
+    hbar_R = hbar * R(np, nh, eom_norm);
   }
 
   // connectivity: empty for unitary ansatz, build otherwise
@@ -361,9 +366,10 @@ std::vector<ExprPtr> CC::eom_r(nₚ np, nₕ nh) const {
   std::int64_t rp = np, rh = nh;
   while (rp >= 0 && rh >= 0) {
     if (rp == 0 && rh == 0) break;
-    // project with <rp, rh| (i.e., multiply P(rp, rh)) and compute VEV
+    // project with <rp, rh| (i.e., multiply δl(rp, rh)) and compute VEV
+    // use δl for consistent normalization
     result.at(min(rp, rh)) =
-        this->ref_av(P(nₚ(rp), nₕ(rh)) * hbar_R, op_connect);
+        this->ref_av(δl(nₚ(rp), nₕ(rh)) * hbar_R, op_connect);
     if (rp == 0 || rh == 0) break;
     rp--;
     rh--;
@@ -386,7 +392,7 @@ std::vector<ExprPtr> CC::eom_l(nₚ np, nₕ nh) const {
   const auto hbar = this->hbar();
 
   // L * hbar
-  const auto L_hbar = L(np, nh) * hbar;
+  const auto L_hbar = L(np, nh, eom_norm) * hbar;
 
   // connectivity:
   // default connections + connect H with projectors
@@ -405,9 +411,9 @@ std::vector<ExprPtr> CC::eom_l(nₚ np, nₕ nh) const {
   std::int64_t rp = np, rh = nh;
   while (rp >= 0 && rh >= 0) {
     if (rp == 0 && rh == 0) break;
-    // right project with |rp,rh> (i.e., multiply P(-rp, -rh)) and compute VEV
+    // right project with |rp,rh> (i.e., multiply δr(rp, rh)) and compute VEV
     result.at(min(rp, rh)) =
-        this->ref_av(L_hbar * P(nₚ(-rp), nₕ(-rh)), op_connect);
+        this->ref_av(L_hbar * δr(nₚ(rp), nₕ(rh)), op_connect);
     if (rp == 0 || rh == 0) break;
     rp--;
     rh--;
diff --git a/SeQuant/domain/mbpt/op.cpp b/SeQuant/domain/mbpt/op.cpp
index c69c20ece5..dc19176dde 100644
--- a/SeQuant/domain/mbpt/op.cpp
+++ b/SeQuant/domain/mbpt/op.cpp
@@ -557,8 +557,9 @@ OpMaker<S>::OpMaker(const std::wstring& label, ncre nc, nann na,
 }
 
 template <Statistics S>
-ExprPtr OpMaker<S>::operator()(std::optional<UseDepIdx> dep,
-                               std::optional<Symmetry> opsymm_opt) const {
+ExprPtr OpMaker<S>::operator()(
+    std::optional<UseDepIdx> dep, std::optional<Symmetry> opsymm_opt,
+    std::optional<Normalization> normalization) const {
   auto isr = get_default_context(Statistics::FermiDirac).index_space_registry();
 
   // if not given dep, use mbpt::Context::CSV to determine whether to use
@@ -597,10 +598,12 @@ ExprPtr OpMaker<S>::operator()(std::optional<UseDepIdx> dep,
   }
   const auto full_label = detail::decorate_with_pert_order(label_, order_);
 
-  const auto normalization =
-      label_ == reserved::antisymm_label() || label_ == reserved::symm_label()
-          ? Normalization::Implicit
-          : Normalization::Default;
+  if (!normalization) {
+    normalization =
+        label_ == reserved::antisymm_label() || label_ == reserved::symm_label()
+            ? Normalization::Implicit
+            : Normalization::Default;
+  }
 
   // if batching indices are present, use them
   if (batch_indices_) {
@@ -613,7 +616,7 @@ ExprPtr OpMaker<S>::operator()(std::optional<UseDepIdx> dep,
                             aux(batchidxs), opsymm_opt ? *opsymm_opt : opsymm,
                             op_herm);
         },
-        dep ? *dep : UseDepIdx::None, normalization);
+        dep ? *dep : UseDepIdx::None, normalization.value());
   }
   // else no batching
   return make(
@@ -623,7 +626,7 @@ ExprPtr OpMaker<S>::operator()(std::optional<UseDepIdx> dep,
         return ex<Tensor>(full_label, bra(creidxs), ket(annidxs),
                           opsymm_opt ? *opsymm_opt : opsymm, op_herm);
       },
-      dep ? *dep : UseDepIdx::None, normalization);
+      dep ? *dep : UseDepIdx::None, normalization.value());
 }
 
 template class OpMaker<Statistics::FermiDirac>;
@@ -756,41 +759,43 @@ ExprPtr Λ(std::size_t K, bool skip1) {
 }
 
 ExprPtr r(nann na, ncre nc, const cre<IndexSpace>& cre_space,
-          const ann<IndexSpace>& ann_space) {
+          const ann<IndexSpace>& ann_space, Normalization norm) {
   SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"R"));
-  return OpMaker<Statistics::FermiDirac>(L"R", nc, na, cre_space, ann_space)();
+  return OpMaker<Statistics::FermiDirac>(L"R", nc, na, cre_space, ann_space)(
+      {}, {}, norm);
 }
-ExprPtr r(nₚ np, nₕ nh) {
+ExprPtr r(nₚ np, nₕ nh, Normalization norm) {
   SEQUANT_ASSERT(np >= 0 && nh >= 0);
   SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"R"));
   return OpMaker<Statistics::FermiDirac>(L"R", ncre(np.value()),
-                                         nann(nh.value()))();
+                                         nann(nh.value()))({}, {}, norm);
 }
 
 ExprPtr l(nann na, ncre nc, const cre<IndexSpace>& cre_space,
-          const ann<IndexSpace>& ann_space) {
+          const ann<IndexSpace>& ann_space, Normalization norm) {
   SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"L"));
-  return OpMaker<Statistics::FermiDirac>(L"L", nc, na, cre_space, ann_space)();
+  return OpMaker<Statistics::FermiDirac>(L"L", nc, na, cre_space, ann_space)(
+      {}, {}, norm);
 }
 
-ExprPtr l(nₚ np, nₕ nh) {
+ExprPtr l(nₚ np, nₕ nh, Normalization norm) {
   SEQUANT_ASSERT(np >= 0 && nh >= 0);
   SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"L"));
   return OpMaker<Statistics::FermiDirac>(L"L", ncre(nh.value()),
-                                         nann(np.value()))();
+                                         nann(np.value()))({}, {}, norm);
 }
 
-ExprPtr P(nₚ np, nₕ nh) {
+ExprPtr P(nₚ np, nₕ nh, std::optional<Normalization> norm) {
   if (np != nh)
     SEQUANT_ASSERT(
         get_default_context().spbasis() != SPBasis::Spinfree &&
         "Spinfree basis does not support non-particle conserving projectors");
   return get_default_context().spbasis() == SPBasis::Spinfree
-             ? tensor::S(-nh /* nh == np */)
-             : tensor::A(-np, -nh);
+             ? tensor::S(-nh /* nh == np */, norm)
+             : tensor::A(-np, -nh, norm);
 }
 
-ExprPtr A(nₚ np, nₕ nh) {
+ExprPtr A(nₚ np, nₕ nh, std::optional<Normalization> norm) {
   SEQUANT_ASSERT(!(np == 0 && nh == 0));
   // if one of them is not zero, nh and np should have the same sign
   if (np != 0 && nh != 0) {
@@ -819,10 +824,10 @@ ExprPtr A(nₚ np, nₕ nh) {
                              : OpMaker<Statistics::FermiDirac>::UseDepIdx::Ket;
   return OpMaker<Statistics::FermiDirac>(reserved::antisymm_label(),
                                          cre(creators), ann(annihilators))(
-      dep, {Symmetry::Antisymm});
+      dep, {Symmetry::Antisymm}, norm);
 }
 
-ExprPtr S(std::int64_t K) {
+ExprPtr S(std::int64_t K, std::optional<Normalization> norm) {
   SEQUANT_ASSERT(K != 0);
   container::svector<IndexSpace> creators;
   container::svector<IndexSpace> annihilators;
@@ -845,7 +850,7 @@ ExprPtr S(std::int64_t K) {
                 : OpMaker<Statistics::FermiDirac>::UseDepIdx::Ket;
   return OpMaker<Statistics::FermiDirac>(reserved::symm_label(), cre(creators),
                                          ann(annihilators))(
-      dep, {Symmetry::Nonsymm});
+      dep, {Symmetry::Nonsymm}, norm);
 }
 
 ExprPtr Hʼ(std::size_t R, const OpParams& params) {
@@ -897,6 +902,19 @@ ExprPtr Λʼ(std::size_t K, const OpParams& params) {
   return result;
 }
 
+// δr/δl are the (de)excitation projectors P, normalized by SquareRoot and
+// indexed by nonnegative ranks (think "derivative wrt r/l"). δl is the
+// deexcitation/bra projector P(np,nh); δr is the excitation/ket projector
+// P(-np,-nh).
+ExprPtr δr(nₚ np, nₕ nh) {
+  SEQUANT_ASSERT(np >= 0 && nh >= 0);
+  return tensor::P(-np, -nh, Normalization::SquareRoot);
+}
+
+ExprPtr δl(nₚ np, nₕ nh) {
+  SEQUANT_ASSERT(np >= 0 && nh >= 0);
+  return tensor::P(np, nh, Normalization::SquareRoot);
+}
 }  // namespace tensor
 
 ExprPtr h(std::size_t k) {
@@ -1013,7 +1031,7 @@ ExprPtr F(bool use_f_tensor, const IndexSpace& occupied_density) {
   }
 }
 
-ExprPtr A(nₚ np, nₕ nh) {
+ExprPtr A(nₚ np, nₕ nh, std::optional<Normalization> norm) {
   SEQUANT_ASSERT(!(nh == 0 && np == 0));
   // if one of them is not zero, nh and np should have the same sign
   if (nh != 0 && np != 0) {
@@ -1026,7 +1044,7 @@ ExprPtr A(nₚ np, nₕ nh) {
   auto hole_space = get_hole_space(Spin::any);
   return ex<op_t>(
       []() -> std::wstring_view { return reserved::antisymm_label(); },
-      [=]() -> ExprPtr { return tensor::A(np, nh); },
+      [=]() -> ExprPtr { return tensor::A(np, nh, norm); },
       [=](qnc_t& qns) {
         const std::size_t abs_nh = std::abs(nh);
         const std::size_t abs_np = std::abs(np);
@@ -1042,10 +1060,10 @@ ExprPtr A(nₚ np, nₕ nh) {
       });
 }
 
-ExprPtr S(std::int64_t K) {
+ExprPtr S(std::int64_t K, std::optional<Normalization> norm) {
   SEQUANT_ASSERT(K != 0);
   return ex<op_t>([]() -> std::wstring_view { return reserved::symm_label(); },
-                  [=]() -> ExprPtr { return tensor::S(K); },
+                  [=]() -> ExprPtr { return tensor::S(K, norm); },
                   [=](qnc_t& qns) {
                     const std::size_t abs_K = std::abs(K);
                     if (K < 0) {
@@ -1058,17 +1076,17 @@ ExprPtr S(std::int64_t K) {
                   });
 }
 
-ExprPtr P(nₚ np, nₕ nh) {
+ExprPtr P(nₚ np, nₕ nh, std::optional<Normalization> norm) {
   if (get_default_context().spbasis() == SPBasis::Spinfree) {
     SEQUANT_ASSERT(
         nh == np &&
         "Only particle number conserving cases are supported with spinfree "
         "basis for now");
     const auto K = np;  // K = np = nh
-    return S(-K);
+    return S(-K, norm);
   } else {
     SEQUANT_ASSERT(get_default_context().spbasis() == SPBasis::Spinor);
-    return A(-np, -nh);
+    return A(-np, -nh, norm);
   }
 }
 
@@ -1129,41 +1147,51 @@ ExprPtr Λʼ(std::size_t K, const OpParams& params) {
 }
 
 ExprPtr r(nann na, ncre nc, const cre<IndexSpace>& cre_space,
-          const ann<IndexSpace>& ann_space) {
+          const ann<IndexSpace>& ann_space, Normalization norm) {
   SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"R"));
-  return ex<op_t>(
-      []() -> std::wstring_view { return L"R"; },
-      [=]() -> ExprPtr { return tensor::r(na, nc, cre_space, ann_space); },
-      [=](qnc_t& qns) {
-        // ex -> creators in particle_space, annihilators in hole_space
-        qns = combine(
-            generic_excitation_qns(/*particle_rank*/ nc, /*hole_rank*/ na,
-                                   cre_space, ann_space),
-            qns);
-      });
+  return ex<op_t>([]() -> std::wstring_view { return L"R"; },
+                  [=]() -> ExprPtr {
+                    return tensor::r(na, nc, cre_space, ann_space, norm);
+                  },
+                  [=](qnc_t& qns) {
+                    // ex -> creators in particle_space, annihilators in
+                    // hole_space
+                    qns = combine(generic_excitation_qns(/*particle_rank*/ nc,
+                                                         /*hole_rank*/ na,
+                                                         cre_space, ann_space),
+                                  qns);
+                  });
 }
 
-ExprPtr r(nₚ np, nₕ nh) { return r(nann(nh), ncre(np)); }
+ExprPtr r(nₚ np, nₕ nh, Normalization norm) {
+  return r(nann(nh), ncre(np), cre(get_particle_space(Spin::any)),
+           ann(get_hole_space(Spin::any)), norm);
+}
 
 ExprPtr l(nann na, ncre nc, const cre<IndexSpace>& cre_space,
-          const ann<IndexSpace>& ann_space) {
+          const ann<IndexSpace>& ann_space, Normalization norm) {
   SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"L"));
-  return ex<op_t>(
-      []() -> std::wstring_view { return L"L"; },
-      [=]() -> ExprPtr { return tensor::l(na, nc, cre_space, ann_space); },
-      [=](qnc_t& qns) {
-        // deex -> creators in hole_space, annihilators in particle_space
-        qns = combine(
-            generic_deexcitation_qns(
-                /*particle_rank*/ na, /*hole_rank*/ nc, ann_space, cre_space),
-            qns);
-      });
+  return ex<op_t>([]() -> std::wstring_view { return L"L"; },
+                  [=]() -> ExprPtr {
+                    return tensor::l(na, nc, cre_space, ann_space, norm);
+                  },
+                  [=](qnc_t& qns) {
+                    // deex -> creators in hole_space, annihilators in
+                    // particle_space
+                    qns = combine(generic_deexcitation_qns(
+                                      /*particle_rank*/ na, /*hole_rank*/ nc,
+                                      ann_space, cre_space),
+                                  qns);
+                  });
 }
 
-ExprPtr l(nₚ np, nₕ nh) { return l(nann(np), ncre(nh)); }
+ExprPtr l(nₚ np, nₕ nh, Normalization norm) {
+  return l(nann(np), ncre(nh), cre(get_hole_space(Spin::any)),
+           ann(get_particle_space(Spin::any)), norm);
+}
 
 ExprPtr R(nann na, ncre nc, const cre<IndexSpace>& cre_space,
-          const ann<IndexSpace>& ann_space) {
+          const ann<IndexSpace>& ann_space, Normalization norm) {
   SEQUANT_ASSERT(na > 0 || nc > 0);
   SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"R"));
   ExprPtr result;
@@ -1171,7 +1199,7 @@ ExprPtr R(nann na, ncre nc, const cre<IndexSpace>& cre_space,
   std::int64_t ra = na, rc = nc;
   while (ra >= 0 && rc >= 0) {
     if (ra == 0 && rc == 0) break;
-    result += r(nann(ra), ncre(rc), cre_space, ann_space);
+    result += r(nann(ra), ncre(rc), cre_space, ann_space, norm);
     if (ra == 0 || rc == 0) break;
     --ra;
     --rc;
@@ -1179,10 +1207,13 @@ ExprPtr R(nann na, ncre nc, const cre<IndexSpace>& cre_space,
   return result;
 }
 
-ExprPtr R(nₚ np, nₕ nh) { return R(nann(nh), ncre(np)); }
+ExprPtr R(nₚ np, nₕ nh, Normalization norm) {
+  return R(nann(nh), ncre(np), cre(get_particle_space(Spin::any)),
+           ann(get_hole_space(Spin::any)), norm);
+}
 
 ExprPtr L(nann na, ncre nc, const cre<IndexSpace>& cre_space,
-          const ann<IndexSpace>& ann_space) {
+          const ann<IndexSpace>& ann_space, Normalization norm) {
   SEQUANT_ASSERT(na > 0 || nc > 0);
   SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"L"));
   ExprPtr result;
@@ -1190,7 +1221,7 @@ ExprPtr L(nann na, ncre nc, const cre<IndexSpace>& cre_space,
   std::int64_t ra = na, rc = nc;
   while (ra >= 0 && rc >= 0) {
     if (ra == 0 && rc == 0) break;
-    result += l(nann(ra), ncre(rc), cre_space, ann_space);
+    result += l(nann(ra), ncre(rc), cre_space, ann_space, norm);
     if (ra == 0 || rc == 0) break;
     --ra;
     --rc;
@@ -1198,7 +1229,24 @@ ExprPtr L(nann na, ncre nc, const cre<IndexSpace>& cre_space,
   return result;
 }
 
-ExprPtr L(nₚ np, nₕ nh) { return L(nann(np), ncre(nh)); }
+ExprPtr L(nₚ np, nₕ nh, Normalization norm) {
+  return L(nann(np), ncre(nh), cre(get_hole_space(Spin::any)),
+           ann(get_particle_space(Spin::any)), norm);
+}
+
+// δr/δl are the (de)excitation projectors P, normalized by SquareRoot and
+// indexed by nonnegative ranks (think "derivative wrt r/l"). δl is the
+// deexcitation/bra projector P(np,nh); δr is the excitation/ket projector
+// P(-np,-nh).
+ExprPtr δr(nₚ np, nₕ nh) {
+  SEQUANT_ASSERT(np >= 0 && nh >= 0);
+  return P(-np, -nh, Normalization::SquareRoot);
+}
+
+ExprPtr δl(nₚ np, nₕ nh) {
+  SEQUANT_ASSERT(np >= 0 && nh >= 0);
+  return P(np, nh, Normalization::SquareRoot);
+}
 
 qns_t apply_to_vac(const ExprPtr& expr) {
   SEQUANT_ASSERT(expr.is<op_t>() || expr.is<Product>());
diff --git a/SeQuant/domain/mbpt/op.hpp b/SeQuant/domain/mbpt/op.hpp
index 993c13baca..da1ce8bfe3 100644
--- a/SeQuant/domain/mbpt/op.hpp
+++ b/SeQuant/domain/mbpt/op.hpp
@@ -567,6 +567,17 @@ mbpt::qns_t adjoint(mbpt::qns_t qns);
 
 namespace mbpt {
 
+/// @brief Normalization convention used in MBPT Operators
+/// An Op of rank {c,a}, where `c`/`a` are the number of creators/annihilators
+/// by default includes a normalization factor of 1/(c! a!).
+/// For some cases, we want to change that.
+enum class Normalization {
+  Default,    /// Include 1/(c! a!) prefactor
+  Implicit,   /// No prefactor, used for  and Ŝ since their definition
+              /// includes the normalization
+  SquareRoot  /// Include sqrt(1/c! a!) prefactor
+};
+
 // clang-format off
 /// @brief makes a tensor-level many-body operator
 
@@ -641,12 +652,6 @@ class OpMaker {
     None
   };
 
-  /// Op of e.g. rank {c,a} by default includes normalization factor
-  /// of 1/(c! a!) ... sometimes we want to change normalization, e.g.
-  /// A and S include normalization factor in their definition,
-  /// so no need to include it explicitly
-  enum class Normalization { Default, Implicit };
-
   /// struct to hold the information about the operator
   struct OpInfo {
     container::svector<Index> creidxs;  //!< creator indices
@@ -666,9 +671,11 @@ class OpMaker {
   /// @param[in] opsymm_opt if given, controls whether (anti)symmetric
   /// tensor is returned; if \p opsymm_opt is not given then the default is
   /// determined by the MBPT context.
+  /// @param[in] normalization if given, controls the normalization behavior, else uses internal defaults. @see Normalization
   // clang-format on
   ExprPtr operator()(std::optional<UseDepIdx> dep_opt = {},
-                     std::optional<Symmetry> opsymm_opt = {}) const;
+                     std::optional<Symmetry> opsymm_opt = {},
+                     std::optional<Normalization> normalization = {}) const;
 
   /// @brief Creates an OpInfo struct containing creator and annihilator
   /// indices, normalization factor, symmetry, and dependency information.
@@ -683,7 +690,7 @@ class OpMaker {
                               const IndexSpaceContainer& ann_spaces,
                               UseDepIdx dep = UseDepIdx::None) {
     const bool symm = get_default_context().spbasis() ==
-                      SPBasis::Spinor;  // antisymmetrize if spin-orbital basis
+                      SPBasis::Spinor;  // antisymmetrize if spinor basis
     const auto dep_bra = dep == UseDepIdx::Bra;
     const auto dep_ket = dep == UseDepIdx::Ket;
 
@@ -729,6 +736,24 @@ class OpMaker {
     return OpInfo{creidxs, annidxs, mult, opsymm, dep};
   }
 
+  /// @brief Applies the prefactor implied by \p normalization (a function of
+  /// the normalization factor \p mult) to \p expr.
+  /// @see Normalization
+  static ExprPtr apply_normalization(ExprPtr expr, Normalization normalization,
+                                     sequant::intmax_t mult) {
+    switch (normalization) {
+      case Normalization::Default:
+        return ex<Constant>(rational{1, mult}) * expr;
+      case Normalization::Implicit:
+        return expr;
+      case Normalization::SquareRoot:
+        return ex<Power>(ex<Constant>(rational{1, mult}), rational{1, 2}) *
+               expr;
+      default:
+        SEQUANT_ABORT("unknown normalization option");
+    }
+  }
+
   /// @tparam TensorGenerator callable with signature
   /// `TensorGenerator(range<Index>, range<Index>, Symmetry)` that returns a
   /// Tensor with the respective bra/cre and ket/ann indices and of the given
@@ -751,16 +776,7 @@ class OpMaker {
     auto result =
         t * ex<NormalOperator<S>>(cre(op_info.creidxs), ann(op_info.annidxs),
                                   get_default_context().vacuum());
-    switch (normalization) {
-      case Normalization::Default:
-        result = ex<Constant>(rational{1, op_info.mult}) * result;
-        break;
-      case Normalization::Implicit:
-        break;
-      default:
-        abort();
-    }
-    return result;
+    return apply_normalization(result, normalization, op_info.mult);
   }
 
   /// @tparam TensorGenerator callable with signature
@@ -818,16 +834,7 @@ class OpMaker {
     auto result =
         t * ex<NormalOperator<S>>(cre(op_info.creidxs), ann(op_info.annidxs),
                                   get_default_context().vacuum());
-    switch (normalization) {
-      case Normalization::Default:
-        result = ex<Constant>(rational{1, op_info.mult}) * result;
-        break;
-      case Normalization::Implicit:
-        break;
-      default:
-        abort();
-    }
-    return result;
+    return apply_normalization(result, normalization, op_info.mult);
   }
 
   /// @tparam TensorGenerator callable with signature
@@ -1014,14 +1021,17 @@ ExprPtr Λ(std::size_t K, bool skip1 = false);
 /// @param nc number of creators
 /// @param cre_space IndexSpace on which creators act
 /// @param ann_space IndexSpace on which annihilators act
+/// @param norm normalization convention, see Normalization
 ExprPtr r(nann na, ncre nc,
           const cre<IndexSpace>& cre_space = cre(get_particle_space(Spin::any)),
-          const ann<IndexSpace>& ann_space = ann(get_hole_space(Spin::any)));
+          const ann<IndexSpace>& ann_space = ann(get_hole_space(Spin::any)),
+          Normalization norm = Normalization::Default);
 
 /// @brief Makes generic excitation operator
 /// @param np number of particle creators
 /// @param nh number of hole creators
-ExprPtr r(nₚ np, nₕ nh);
+/// @param norm normalization convention, see Normalization
+ExprPtr r(nₚ np, nₕ nh, Normalization norm = Normalization::Default);
 DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(r);
 
 /// @brief Makes generic left-hand replacement operator
@@ -1029,25 +1039,28 @@ DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(r);
 /// @param nc number of creators
 /// @param cre_space IndexSpace on which creators act
 /// @param ann_space IndexSpace on which annihilators act
-ExprPtr l(
-    nann na, ncre nc,
-    const cre<IndexSpace>& cre_space = cre(get_hole_space(Spin::any)),
-    const ann<IndexSpace>& ann_space = ann(get_particle_space(Spin::any)));
+/// @param norm normalization convention, see Normalization
+ExprPtr l(nann na, ncre nc,
+          const cre<IndexSpace>& cre_space = cre(get_hole_space(Spin::any)),
+          const ann<IndexSpace>& ann_space = ann(get_particle_space(Spin::any)),
+          Normalization norm = Normalization::Default);
 
 /// @brief Makes generic deexcitation operator
 /// @param np number of particle annihilators
 /// @param nh number of hole annihilators
-ExprPtr l(nₚ np, nₕ nh);
+/// @param norm normalization convention, see Normalization
+ExprPtr l(nₚ np, nₕ nh, Normalization norm = Normalization::Default);
 DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(l);
 
 // clang-format off
 /// makes projector onto excited bra (if \p np > 0 && \p nh > 0) or ket (if \p np < 0 && \p nh <0) manifold
 /// @param np number of particle creators (if > 0) or annihilators (< 0)
 /// @param nh number of hole creators (if > 0) or annihilators (< 0); if omitted, will use \p np
+/// @param norm normalization convention; if unset, uses the intrinsic Implicit normalization. @see Normalization
 /// @note if using spin-free basis, only supports particle-symmetric operators `K = Kh = Kp`, returns `S(-K)`
 /// else supports particle non-conserving operators and returns `A(-np, -nh)`
 // clang-format on
-ExprPtr P(nₚ np, nₕ nh);
+ExprPtr P(nₚ np, nₕ nh, std::optional<Normalization> norm = {});
 DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(P);
 
 // clang-format off
@@ -1055,14 +1068,17 @@ DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(P);
 /// @param np number of particle creators (if > 0) or annihilators (< 0)
 /// @param nh number of hole creators (if > 0) or annihilators (< 0); if omitted, will use \p np
 /// (default is to set \p np to \p nh)
+/// @param norm normalization convention; if unset, uses the intrinsic Implicit normalization. @see Normalization
 /// @note supports particle non-conserving operators
 // clang-format on
-ExprPtr A(nₚ np, nₕ nh);
+ExprPtr A(nₚ np, nₕ nh, std::optional<Normalization> norm = {});
 DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(A);
 
 /// @brief makes generic particle-symmetric excitation (if \p K > 0) or
 /// deexcitation (if \p K < 0) operator of rank `|K|`
-ExprPtr S(std::int64_t K);
+/// @param norm normalization convention; if unset, uses the intrinsic Implicit
+/// normalization. @see Normalization
+ExprPtr S(std::int64_t K, std::optional<Normalization> norm = {});
 
 /// @brief Makes perturbation operator
 /// @param R rank of the perturbation operator
@@ -1102,6 +1118,27 @@ ExprPtr λʼ(std::size_t K, const OpParams& params = {.order = 1});
 /// @pre If batching is used, ISR must contain batching space
 ExprPtr Λʼ(std::size_t K,
            const OpParams& params = {.order = 1, .skip1 = false});
+
+// clang-format off
+/// @brief Makes projector 1/√(np! nh!) A_{a1 a2 ... a_np}^{i1 i2 ... i_nh} a_{i1 i2 ... i_nh}^{a1 a2 ... a_np} (excitation operator).
+/// Unlike P, uses SquareRoot normalization (includes the 1/√(np! nh!) prefactor).
+/// @param np number of particle creators
+/// @param nh number of hole annihilators
+/// @note if using spin-free basis, only supports particle-number-conserving operators (\p np == \p nh)
+// clang-format on
+ExprPtr δr(nₚ np, nₕ nh);
+DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(δr);
+
+// clang-format off
+/// @brief Makes projector 1/√(np! nh!) A^{a1 a2 ... a_np}_{i1 i2 ... i_nh} a^{i1 i2 ... i_nh}_{a1 a2 ... a_np} (deexcitation operator).
+/// Unlike P, uses SquareRoot normalization (includes the 1/√(np! nh!) prefactor).
+/// @param np number of particle annihilators
+/// @param nh number of hole creators
+/// @note if using spin-free basis, only supports particle-number-conserving operators (\p np == \p nh)
+// clang-format on
+ExprPtr δl(nₚ np, nₕ nh);
+DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(δl);
+
 }  // namespace tensor
 }  // namespace op
 
@@ -1146,14 +1183,17 @@ ExprPtr Λ(std::size_t K, bool skip1 = false);
 /// @param nc number of creators
 /// @param cre_space IndexSpace on which creators act
 /// @param ann_space IndexSpace on which annihilators act
+/// @param norm normalization convention, see Normalization
 ExprPtr r(nann na, ncre nc,
           const cre<IndexSpace>& cre_space = cre(get_particle_space(Spin::any)),
-          const ann<IndexSpace>& ann_space = ann(get_hole_space(Spin::any)));
+          const ann<IndexSpace>& ann_space = ann(get_hole_space(Spin::any)),
+          Normalization norm = Normalization::Default);
 
 /// @brief Makes generic excitation operator
 /// @param np number of particle creators
 /// @param nh number of hole creators
-ExprPtr r(nₚ np, nₕ nh);
+/// @param norm normalization convention, see Normalization
+ExprPtr r(nₚ np, nₕ nh, Normalization norm = Normalization::Default);
 DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(r);
 
 /// @brief Makes generic deexcitation operator
@@ -1161,15 +1201,17 @@ DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(r);
 /// @param nc number of creators
 /// @param cre_space IndexSpace on which creators act
 /// @param ann_space IndexSpace on which annihilators act
-ExprPtr l(
-    nann na, ncre nc,
-    const cre<IndexSpace>& cre_space = cre(get_hole_space(Spin::any)),
-    const ann<IndexSpace>& ann_space = ann(get_particle_space(Spin::any)));
+/// @param norm normalization convention, see Normalization
+ExprPtr l(nann na, ncre nc,
+          const cre<IndexSpace>& cre_space = cre(get_hole_space(Spin::any)),
+          const ann<IndexSpace>& ann_space = ann(get_particle_space(Spin::any)),
+          Normalization norm = Normalization::Default);
 
 /// @brief Makes generic deexcitation operator
 /// @param np number of particle annihilators
 /// @param nh number of hole annihilators
-ExprPtr l(nₚ np, nₕ nh);
+/// @param norm normalization convention, see Normalization
+ExprPtr l(nₚ np, nₕ nh, Normalization norm = Normalization::Default);
 DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(l);
 
 /// @brief Makes sum of generic right-hand replacement operators up to max rank
@@ -1177,16 +1219,19 @@ DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(l);
 /// @param nc number of creators
 /// @param cre_space IndexSpace on which creators act
 /// @param ann_space IndexSpace on which annihilators act
+/// @param norm normalization convention, see Normalization
 /// @return `r(na,nc) + r(na-1,nc-1) + ...`
 ExprPtr R(nann na, ncre nc,
           const cre<IndexSpace>& cre_space = cre(get_particle_space(Spin::any)),
-          const ann<IndexSpace>& ann_space = ann(get_hole_space(Spin::any)));
+          const ann<IndexSpace>& ann_space = ann(get_hole_space(Spin::any)),
+          Normalization norm = Normalization::Default);
 
 /// @brief Makes sum of generic excitation operators up to max rank
 /// @param np max number of particle creators
 /// @param nh max number of hole creators
+/// @param norm normalization convention, see Normalization
 /// @return `r(np,nh) + r(np-1,nh-1) + ...`
-ExprPtr R(nₚ np, nₕ nh);
+ExprPtr R(nₚ np, nₕ nh, Normalization norm = Normalization::Default);
 DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(R);
 
 /// @brief Makes sum of generic "left-hand" replacement operators up to max rank
@@ -1194,41 +1239,47 @@ DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(R);
 /// @param nc number of creators
 /// @param cre_space IndexSpace on which creators act
 /// @param ann_space IndexSpace on which annihilators act
+/// @param norm normalization convention, see Normalization
 /// @return `l(na,nc) + l(na-1,nc-1) + ...`
-ExprPtr L(
-    nann na, ncre nc,
-    const cre<IndexSpace>& cre_space = cre(get_hole_space(Spin::any)),
-    const ann<IndexSpace>& ann_space = ann(get_particle_space(Spin::any)));
+ExprPtr L(nann na, ncre nc,
+          const cre<IndexSpace>& cre_space = cre(get_hole_space(Spin::any)),
+          const ann<IndexSpace>& ann_space = ann(get_particle_space(Spin::any)),
+          Normalization norm = Normalization::Default);
 
 /// @brief Makes sum of deexcitation operators up to max rank
 /// @param np max number of particle annihilators
 /// @param nh max number of hole annihilators
+/// @param norm normalization convention, see Normalization
 /// @return `l(np,nh) + l(np-1,nh-1) + ...`
-ExprPtr L(nₚ np, nₕ nh);
+ExprPtr L(nₚ np, nₕ nh, Normalization norm = Normalization::Default);
 DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(L);
 
 // clang-format off
 /// makes projector onto excited bra (if \p np > 0 && \p nh > 0) or ket (if \p np < 0 && \p nh <0) manifold
 /// @param np number of particle creators (if > 0) or annihilators (< 0)
 /// @param nh number of hole creators (if > 0) or annihilators (< 0); if omitted, will use \p np
+/// @param norm normalization convention; if unset, uses the intrinsic Implicit normalization. @see Normalization
 /// @note if using spin-free basis, only supports particle-symmetric operators `K = Kh = Kp`, returns `S(-K)`
 /// else supports particle non-conserving operators and returns `A(-np, -nh)`
 // clang-format on
-ExprPtr P(nₚ np, nₕ nh);
+ExprPtr P(nₚ np, nₕ nh, std::optional<Normalization> norm = {});
 DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(P);
 
 // clang-format off
 /// @brief makes generic bra/ket-antisymmetric excitation (if \p nh > 0 && \p np > 0) or deexcitation (if \p nh < 0 && \p np < 0) operator
 /// @param np number of particle creators (if > 0) or annihilators (< 0)
 /// @param nh number of hole creators (if > 0) or annihilators (< 0); if omitted, will use \p np
+/// @param norm normalization convention; if unset, uses the intrinsic Implicit normalization. @see Normalization
 /// @note supports particle non-conserving operators
 // clang-format on
-ExprPtr A(nₚ np, nₕ nh);
+ExprPtr A(nₚ np, nₕ nh, std::optional<Normalization> norm = {});
 DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(A);
 
 /// @brief makes generic particle-symmetric excitation (if \p K > 0) or
 /// deexcitation (if \p K < 0) operator of rank `|K|`
-ExprPtr S(std::int64_t K);
+/// @param norm normalization convention; if unset, uses the intrinsic Implicit
+/// normalization. @see Normalization
+ExprPtr S(std::int64_t K, std::optional<Normalization> norm = {});
 
 /// @brief Makes perturbation operator
 /// @param R rank of the perturbation operator
@@ -1269,6 +1320,26 @@ ExprPtr λʼ(std::size_t K, const OpParams& params = {.order = 1});
 ExprPtr Λʼ(std::size_t K,
            const OpParams& params = {.order = 1, .skip1 = false});
 
+// clang-format off
+/// @brief Makes projector with tensor form 1/√(np! nh!) A_{a1 a2 ... a_np}^{i1 i2 ... i_nh} a_{i1 i2 ... i_nh}^{a1 a2 ... a_np} (excitation operator).
+/// Unlike P, uses SquareRoot normalization (includes the 1/√(np! nh!) prefactor).
+/// @param np number of particle creators
+/// @param nh number of hole annihilators
+/// @note if using spin-free basis, only supports particle-number-conserving operators (\p np == \p nh)
+// clang-format on
+ExprPtr δr(nₚ np, nₕ nh);
+DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(δr);
+
+// clang-format off
+/// @brief Makes projector with tensor form 1/√(np! nh!) A^{a1 a2 ... a_np}_{i1 i2 ... i_nh} a^{i1 i2 ... i_nh}_{a1 a2 ... a_np} (deexcitation operator).
+/// Unlike P, uses SquareRoot normalization (includes the 1/√(np! nh!) prefactor).
+/// @param np number of particle annihilators
+/// @param nh number of hole creators
+/// @note if using spin-free basis, only supports particle-number-conserving operators (\p np == \p nh)
+// clang-format on
+ExprPtr δl(nₚ np, nₕ nh);
+DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(δl);
+
 /// @brief computes the quantum number change effected by a given Operator or
 /// Operator Product when applied to the vacuum state
 /// @param expr the operator or operator product whose quantum number change is
diff --git a/tests/unit/test_mbpt.cpp b/tests/unit/test_mbpt.cpp
index 0e3e8f4e10..92d7677dc4 100644
--- a/tests/unit/test_mbpt.cpp
+++ b/tests/unit/test_mbpt.cpp
@@ -781,6 +781,34 @@ TEST_CASE("mbpt", "[mbpt][valgrind_skip]") {
 
       REQUIRE(to_latex(simplify(h1 * t2)) == L"{{\\hat{h¹}}{\\hat{t²}_{2}}}");
       REQUIRE(to_latex(simplify(h2 * t2)) == L"{{\\hat{h²}}{\\hat{t²}_{2}}}");
+
+      // δl δr ops
+      {  // Spinor basis
+        auto dl2 = tensor::δl(2);
+        REQUIRE(simplify(dl2 - rational{1, 2} * tensor::P(2)) ==
+                ex<Constant>(0));
+        auto dr2 = tensor::δr(2);
+        REQUIRE(simplify(dr2 - rational{1, 2} * tensor::P(-2)) ==
+                ex<Constant>(0));
+
+        auto dl23 = tensor::δl(nₚ(2), nₕ(3));
+        REQUIRE(simplify(dl23 - ex<Power>(rational{1, 12}, rational{1, 2}) *
+                                    tensor::P(nₚ(2), nₕ(3))) ==
+                ex<Constant>(0));
+      }
+
+      {  // spinfree basis
+        auto ctx = get_default_context();
+        auto ctx_resetter =
+            set_scoped_default_context(ctx.set(SPBasis::Spinfree));
+        auto dl2 = tensor::δl(2);
+        REQUIRE(simplify(dl2 - ex<Power>(rational{1, 2}, rational{1, 2}) *
+                                   tensor::P(2)) == ex<Constant>(0));
+        auto dr3 = tensor::δr(3);
+        REQUIRE(simplify(dr3 - ex<Power>(rational{1, 6}, rational{1, 2}) *
+                                   tensor::P(-3)) == ex<Constant>(0));
+      }
+
     }  // SECTION("predefined")
 
     SECTION("batching") {