diff --git a/examples/fourier_transform_1d.cpp b/examples/fourier_transform_1d.cpp
index cf758ef..7961e9f 100644
--- a/examples/fourier_transform_1d.cpp
+++ b/examples/fourier_transform_1d.cpp
@@ -24,8 +24,8 @@ class TightBindingModel1D {
   Expression m_hamiltonian;
 
   void construct_hamiltonian() {
-    for (uint8_t i = 0; i < m_size; ++i) {
-      uint8_t j = (i + 1) % m_size;  // Periodic boundary conditions
+    for (size_t i = 0; i < m_size; ++i) {
+      size_t j = (i + 1) % m_size;  // Periodic boundary conditions
       m_hamiltonian += Expression::hopping(i, j, Operator::Spin::Up) * m_t;
       m_hamiltonian += Expression::hopping(i, j, Operator::Spin::Down) * m_t;
     }
@@ -44,7 +44,7 @@ static Expression fourier_transform_operator(const Operator& op, size_t size) {
         0.0f, -type_sign * factor * static_cast<float>(k * op.orbital()));
     coefficient = std::exp(coefficient) / std::sqrt(static_cast<float>(size));
 
-    Operator transformed_op(op.type(), op.spin(), static_cast<uint8_t>(k));
+    Operator transformed_op(op.type(), op.spin(), k);
     result += Term(coefficient, {transformed_op});
   }
 
diff --git a/examples/fourier_transform_2d.cpp b/examples/fourier_transform_2d.cpp
index cdb334c..1550f08 100644
--- a/examples/fourier_transform_2d.cpp
+++ b/examples/fourier_transform_2d.cpp
@@ -17,14 +17,14 @@ class Index2D {
     }
   }
 
-  uint8_t to_orbital(size_t x, size_t y) const {
+  size_t to_orbital(size_t x, size_t y) const {
     if (x >= m_width || y >= m_height) {
       throw std::out_of_range("Coordinates out of lattice bounds");
     }
-    return static_cast<uint8_t>(y * m_width + x);
+    return y * m_width + x;
   }
 
-  std::pair<size_t, size_t> from_orbital(uint8_t orbital) const {
+  std::pair<size_t, size_t> from_orbital(size_t orbital) const {
     if (orbital >= m_width * m_height) {
       throw std::out_of_range("Orbital index out of lattice bounds");
     }
@@ -55,10 +55,10 @@ class TightBindingModel2D {
   Expression m_hamiltonian;
 
   void construct_hamiltonian() {
-    for (uint8_t y = 0; y < m_index.height(); ++y) {
-      for (uint8_t x = 0; x < m_index.width(); ++x) {
+    for (size_t y = 0; y < m_index.height(); ++y) {
+      for (size_t x = 0; x < m_index.width(); ++x) {
         // Horizontal hopping
-        uint8_t next_x = (x + 1) % m_index.width();
+        size_t next_x = (x + 1) % m_index.width();
         m_hamiltonian +=
             m_t * Expression::hopping(m_index.to_orbital(x, y),
                                       m_index.to_orbital(next_x, y),
@@ -69,7 +69,7 @@ class TightBindingModel2D {
                                       Operator::Spin::Down);
 
         // Vertical hopping
-        uint8_t next_y = (y + 1) % m_index.height();
+        size_t next_y = (y + 1) % m_index.height();
         m_hamiltonian +=
             m_t * Expression::hopping(m_index.to_orbital(x, y),
                                       m_index.to_orbital(x, next_y),
diff --git a/examples/hubbard_1d.cpp b/examples/hubbard_1d.cpp
index ca202d3..b469280 100644
--- a/examples/hubbard_1d.cpp
+++ b/examples/hubbard_1d.cpp
@@ -27,8 +27,8 @@ class HubbardChain1D {
 
 void HubbardChain1D::construct_hamiltonian() {
   // Hopping terms
-  for (uint8_t i = 0; i < m_sites; ++i) {
-    uint8_t j = (i + 1) % m_sites;  // Periodic boundary conditions
+  for (size_t i = 0; i < m_sites; ++i) {
+    size_t j = (i + 1) % m_sites;  // Periodic boundary conditions
     m_hamiltonian +=
         qmutils::Expression::hopping(i, j, qmutils::Operator::Spin::Up);
     m_hamiltonian +=
@@ -36,7 +36,7 @@ void HubbardChain1D::construct_hamiltonian() {
   }
 
   // On-site interaction terms
-  for (uint8_t i = 0; i < m_sites; ++i) {
+  for (size_t i = 0; i < m_sites; ++i) {
     m_hamiltonian += qmutils::Term::density(qmutils::Operator::Spin::Up, i) *
                      qmutils::Term::density(qmutils::Operator::Spin::Down, i) *
                      m_U;
diff --git a/examples/hubbard_1d_mean_field.cpp b/examples/hubbard_1d_mean_field.cpp
index c701ea7..d33821d 100644
--- a/examples/hubbard_1d_mean_field.cpp
+++ b/examples/hubbard_1d_mean_field.cpp
@@ -20,8 +20,8 @@ class Hubbard1DModel {
     qmutils::Expression H;
 
     // Hopping terms
-    for (uint8_t i = 0; i < m_sites; ++i) {
-      uint8_t j = (i + 1) % m_sites;  // Periodic boundary conditions
+    for (size_t i = 0; i < m_sites; ++i) {
+      size_t j = (i + 1) % m_sites;  // Periodic boundary conditions
       H += -m_t *
            qmutils::Expression::hopping(i, j, qmutils::Operator::Spin::Up);
       H += -m_t *
@@ -29,7 +29,7 @@ class Hubbard1DModel {
     }
 
     // Mean-field interaction terms
-    for (uint8_t i = 0; i < m_sites; ++i) {
+    for (size_t i = 0; i < m_sites; ++i) {
       H += m_U * down_occupations[i] *
            qmutils::Term::density(qmutils::Operator::Spin::Up, i);
       H += m_U * up_occupations[i] *
diff --git a/examples/ssh_model.cpp b/examples/ssh_model.cpp
index 67263e2..6b52cdb 100644
--- a/examples/ssh_model.cpp
+++ b/examples/ssh_model.cpp
@@ -25,8 +25,8 @@ class SSHModel {
     auto hopping = [this](size_t unit_cell_1, size_t size_1, size_t unit_cell_2,
                           size_t site_2, float t) {
       Expression term;
-      uint8_t orbital1 = m_index.to_orbital(unit_cell_1, size_1);
-      uint8_t orbital2 = m_index.to_orbital(unit_cell_2, site_2);
+      size_t orbital1 = m_index.to_orbital(unit_cell_1, size_1);
+      size_t orbital2 = m_index.to_orbital(unit_cell_2, site_2);
       term += t * Expression::hopping(orbital1, orbital2, Operator::Spin::Up);
       term += t * Expression::hopping(orbital1, orbital2, Operator::Spin::Down);
       return term;
@@ -35,12 +35,12 @@ class SSHModel {
     const size_t L = m_index.dimension(0);
 
     // Intra-cell hopping
-    for (uint8_t i = 0; i < L; ++i) {
+    for (size_t i = 0; i < L; ++i) {
       H += hopping(i, 0, i, 1, m_t1 + m_delta);
     }
 
     // Inter-cell hopping
-    for (uint8_t i = 0; i < L; ++i) {
+    for (size_t i = 0; i < L; ++i) {
       H += hopping(i, 1, (i + 1) % L, 0, m_t2 - m_delta);
     }
 
@@ -100,7 +100,7 @@ static Expression transform_operator_to_band_basis(
                                     : std::conj(eigenvectors(i, k));
     // Get spin and orbital indices from basis index k
     Operator::Spin spin = basis.at(k)[0].spin();
-    uint8_t orbital = basis.at(k)[0].orbital();
+    size_t orbital = basis.at(k)[0].orbital();
     Term new_term(coeff, {Operator(op.type(), spin, orbital)});
     result += new_term;
   }
diff --git a/src/expression.cpp b/src/expression.cpp
index c0c0885..ef077c6 100644
--- a/src/expression.cpp
+++ b/src/expression.cpp
@@ -62,7 +62,7 @@ Expression Expression::flip_spin() const {
   return result;
 }
 
-Expression Expression::hopping(uint8_t from_orbital, uint8_t to_orbital,
+Expression Expression::hopping(size_t from_orbital, size_t to_orbital,
                                Operator::Spin spin) {
   QMUTILS_ASSERT(from_orbital != to_orbital);
   Expression result;
diff --git a/src/qmutils/expression.h b/src/qmutils/expression.h
index a6e73a6..c73de0b 100644
--- a/src/qmutils/expression.h
+++ b/src/qmutils/expression.h
@@ -167,7 +167,7 @@ class Expression {
 
   Expression flip_spin() const;
 
-  static Expression hopping(uint8_t from_orbital, uint8_t to_orbital,
+  static Expression hopping(size_t from_orbital, size_t to_orbital,
                             Operator::Spin spin);
 
   struct Spin;
@@ -177,7 +177,7 @@ class Expression {
 };
 
 struct Expression::Spin {
-  static Expression spin_x(uint8_t i) {
+  static Expression spin_x(size_t i) {
     Expression result;
     result +=
         0.5f * Term::one_body(Operator::Spin::Up, i, Operator::Spin::Down, i);
@@ -186,7 +186,7 @@ struct Expression::Spin {
     return result;
   }
 
-  static Expression spin_y(uint8_t i) {
+  static Expression spin_y(size_t i) {
     Expression result;
     result += std::complex<float>(0.0f, 0.5f) *
               Term::one_body(Operator::Spin::Up, i, Operator::Spin::Down, i);
@@ -195,7 +195,7 @@ struct Expression::Spin {
     return result;
   }
 
-  static Expression spin_z(uint8_t i) {
+  static Expression spin_z(size_t i) {
     Expression result;
     result +=
         0.5f * Term::one_body(Operator::Spin::Up, i, Operator::Spin::Up, i);
@@ -204,17 +204,17 @@ struct Expression::Spin {
     return result;
   }
 
-  static Expression spin_plus(uint8_t site) {
+  static Expression spin_plus(size_t site) {
     return Expression(
         Term::one_body(Operator::Spin::Up, site, Operator::Spin::Down, site));
   }
 
-  static Expression spin_minus(uint8_t site) {
+  static Expression spin_minus(size_t site) {
     return Expression(
         Term::one_body(Operator::Spin::Down, site, Operator::Spin::Up, site));
   }
 
-  static Expression dot_product(uint8_t i, uint8_t j) {
+  static Expression dot_product(size_t i, size_t j) {
     Expression result;
     result += spin_x(i) * spin_x(j);
     result += spin_y(i) * spin_y(j);
diff --git a/src/qmutils/operator.h b/src/qmutils/operator.h
index 78f533c..a0f22ff 100644
--- a/src/qmutils/operator.h
+++ b/src/qmutils/operator.h
@@ -5,6 +5,8 @@
 #include <string>
 #include <type_traits>
 
+#include "assert.h"
+
 namespace qmutils {
 
 class Operator {
@@ -28,7 +30,9 @@ class Operator {
       : data_{.orbital = static_cast<int_type>(orbital),
               .spin = static_cast<int_type>(spin),
               .type = static_cast<int_type>(type),
-              .statistics = static_cast<int_type>(stat)} {}
+              .statistics = static_cast<int_type>(stat)} {
+    QMUTILS_ASSERT(orbital <= max_orbital_size());
+  }
 
   [[nodiscard]] constexpr Type type() const noexcept {
     return static_cast<Type>(data_.type);
diff --git a/src/qmutils/term.h b/src/qmutils/term.h
index 8c6f123..3a086a9 100644
--- a/src/qmutils/term.h
+++ b/src/qmutils/term.h
@@ -83,11 +83,11 @@ class Term {
 
   std::string to_string() const;
 
-  static Term creation(Operator::Spin spin, uint8_t orbital);
-  static Term annihilation(Operator::Spin spin, uint8_t orbital);
-  static Term one_body(Operator::Spin spin1, uint8_t orbital1,
-                       Operator::Spin spin2, uint8_t orbital2);
-  static Term density(Operator::Spin spin, uint8_t orbital);
+  static Term creation(Operator::Spin spin, size_t orbital);
+  static Term annihilation(Operator::Spin spin, size_t orbital);
+  static Term one_body(Operator::Spin spin1, size_t orbital1,
+                       Operator::Spin spin2, size_t orbital2);
+  static Term density(Operator::Spin spin, size_t orbital);
 
  private:
   coefficient_type m_coefficient;
diff --git a/src/term.cpp b/src/term.cpp
index 499ca2e..62c17f2 100644
--- a/src/term.cpp
+++ b/src/term.cpp
@@ -34,21 +34,21 @@ std::string Term::to_string() const {
   return oss.str();
 }
 
-Term Term::creation(Operator::Spin spin, uint8_t orbital) {
+Term Term::creation(Operator::Spin spin, size_t orbital) {
   return Term(1.0f, {Operator::creation(spin, orbital)});
 }
 
-Term Term::annihilation(Operator::Spin spin, uint8_t orbital) {
+Term Term::annihilation(Operator::Spin spin, size_t orbital) {
   return Term(1.0f, {Operator::annihilation(spin, orbital)});
 }
 
-Term Term::one_body(Operator::Spin spin1, uint8_t orbital1,
-                    Operator::Spin spin2, uint8_t orbital2) {
+Term Term::one_body(Operator::Spin spin1, size_t orbital1, Operator::Spin spin2,
+                    size_t orbital2) {
   return Term(1.0f, {Operator::creation(spin1, orbital1),
                      Operator::annihilation(spin2, orbital2)});
 }
 
-Term Term::density(Operator::Spin spin, uint8_t orbital) {
+Term Term::density(Operator::Spin spin, size_t orbital) {
   return Term::one_body(spin, orbital, spin, orbital);
 }