diff --git a/CMakeLists.txt b/CMakeLists.txt
index ee965c40..c1fd2214 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -32,8 +32,12 @@ include(nwx_cxx_api_docs)
 nwx_cxx_api_docs("cxx/include" "cxx/src")
 
 # Dependencies
+set(TENSORWRAPPER_DEPENDENCIES utilities parallelzone wtf eigen boost)
+if(ENABLE_SIGMA)
+    list(APPEND TENSORWRAPPER_DEPENDENCIES sigma)
+endif()
 include(get_dependencies)
-get_dependencies(utilities parallelzone wtf eigen boost)
+get_dependencies(${TENSORWRAPPER_DEPENDENCIES})
 
 # TODO: ENABLE_CUTENSOR handling
 
@@ -43,6 +47,26 @@ nwx_library(${PROJECT_NAME} "cxx/include" "cxx/src"
     utilities parallelzone wtf
     PUBLIC eigen boost)
 
+# ENABLE_SIGMA is declared as an option() by set_default_nwx_options, but that
+# alone doesn't reach the preprocessor -- floating_point.hpp's #ifdef
+# ENABLE_SIGMA needs the compile definition too.
+target_compile_definitions(${PROJECT_NAME} PUBLIC
+    $<$<BOOL:${ENABLE_SIGMA}>:ENABLE_SIGMA>)
+
+# sigma is linked BUILD_INTERFACE-only: it's still on its own CMaize-based
+# build (see NWXCMake/cmake/dependencies/sigma.cmake), so the raw "sigma"
+# target FetchContent produces is not IMPORTED and isn't part of this
+# project's own install(EXPORT) set. Requiring it there fails at generate
+# time ("target tensorwrapper... requires target sigma that is not in any
+# export set"). BUILD_INTERFACE keeps it available to this build (this
+# target's own sources, and other in-build consumers like the test binaries
+# and Python module) without making the installed/exported tensorwrapper
+# package require it.
+if(ENABLE_SIGMA)
+    target_link_libraries(${PROJECT_NAME} PUBLIC
+        $<BUILD_INTERFACE:${NWX_DEP_TARGET_sigma}>)
+endif()
+
 # Generated config headers (e.g. from configure_file) live in binary dir
 target_include_directories(${PROJECT_NAME} PUBLIC
     $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>)
diff --git a/cxx/include/tensorwrapper/types/floating_point.hpp b/cxx/include/tensorwrapper/types/floating_point.hpp
index dd34104f..fbd75c14 100644
--- a/cxx/include/tensorwrapper/types/floating_point.hpp
+++ b/cxx/include/tensorwrapper/types/floating_point.hpp
@@ -45,9 +45,14 @@ using thresholded_affine_type = sigma::ThresholdedAffine<T>;
 using tafloat                 = thresholded_affine_type<float>;
 using tadouble                = thresholded_affine_type<double>;
 
+template<typename T>
+using taylor_model_type = sigma::TaylorModel<T>;
+using tmfloat           = taylor_model_type<float>;
+using tmdouble          = taylor_model_type<double>;
+
 using floating_point_types =
   std::tuple<float, double, ufloat, udouble, ifloat, idouble, afloat, adouble,
-             tafloat, tadouble>;
+             tafloat, tadouble, tmfloat, tmdouble>;
 
 template<typename T>
 constexpr bool is_uncertain_v =
@@ -66,8 +71,13 @@ constexpr bool is_thresholded_affine_v =
   std::is_same_v<T, tafloat> || std::is_same_v<T, tadouble>;
 
 template<typename T>
-constexpr bool is_uq_type_v = is_uncertain_v<T> || is_interval_v<T> ||
-                              is_affine_v<T> || is_thresholded_affine_v<T>;
+constexpr bool is_taylor_model_v =
+  std::is_same_v<T, tmfloat> || std::is_same_v<T, tmdouble>;
+
+template<typename T>
+constexpr bool is_uq_type_v =
+  is_uncertain_v<T> || is_interval_v<T> || is_affine_v<T> ||
+  is_thresholded_affine_v<T> || is_taylor_model_v<T>;
 
 template<typename T, typename U, typename V>
 T construct_uq_type(const U& center, const V& radius) {
@@ -75,7 +85,8 @@ T construct_uq_type(const U& center, const V& radius) {
         return T(center, radius);
     } else if constexpr(is_interval_v<T>) {
         return T(center - radius, center + radius);
-    } else if constexpr(is_affine_v<T> || is_thresholded_affine_v<T>) {
+    } else if constexpr(is_affine_v<T> || is_thresholded_affine_v<T> ||
+                        is_taylor_model_v<T>) {
         return T(center - radius, center + radius);
     } else if constexpr(is_uq_type_v<T>) {
         throw std::logic_error("UQ type not recognized in construct_uq_type.");
@@ -92,6 +103,8 @@ auto uq_center(const T& value) {
         return value.median();
     } else if constexpr(is_affine_v<T> || is_thresholded_affine_v<T>) {
         return value.center();
+    } else if constexpr(is_taylor_model_v<T>) {
+        return value.constant();
     } else if constexpr(is_uq_type_v<T>) {
         throw std::logic_error("UQ type not recognized in uq_center.");
     } else {
@@ -105,7 +118,8 @@ auto uq_upper(const T& value) {
         return value.mean() + value.sd();
     } else if constexpr(is_interval_v<T>) {
         return value.upper();
-    } else if constexpr(is_affine_v<T> || is_thresholded_affine_v<T>) {
+    } else if constexpr(is_affine_v<T> || is_thresholded_affine_v<T> ||
+                        is_taylor_model_v<T>) {
         return value.range().upper();
     } else if constexpr(is_uq_type_v<T>) {
         throw std::logic_error("UQ type not recognized in uq_upper.");
@@ -114,6 +128,22 @@ auto uq_upper(const T& value) {
     }
 }
 
+template<typename T>
+auto uq_lower(const T& value) {
+    if constexpr(is_uncertain_v<T>) {
+        return value.mean() - value.sd();
+    } else if constexpr(is_interval_v<T>) {
+        return value.lower();
+    } else if constexpr(is_affine_v<T> || is_thresholded_affine_v<T> ||
+                        is_taylor_model_v<T>) {
+        return value.range().lower();
+    } else if constexpr(is_uq_type_v<T>) {
+        throw std::logic_error("UQ type not recognized in uq_lower.");
+    } else {
+        return value;
+    }
+}
+
 template<typename T, typename U>
 bool strictly_less(const T& lhs, const U& rhs) {
     return uq_upper(lhs) < uq_upper(rhs);
@@ -165,7 +195,9 @@ T pow(T value, double pow) {
     MACRO_IN(tensorwrapper::types::afloat);     \
     MACRO_IN(tensorwrapper::types::adouble);    \
     MACRO_IN(tensorwrapper::types::tafloat);    \
-    MACRO_IN(tensorwrapper::types::tadouble)
+    MACRO_IN(tensorwrapper::types::tadouble);   \
+    MACRO_IN(tensorwrapper::types::tmfloat);    \
+    MACRO_IN(tensorwrapper::types::tmdouble)
 } // namespace tensorwrapper::types
 
 WTF_REGISTER_FP_TYPE(tensorwrapper::types::ufloat);
@@ -176,6 +208,8 @@ WTF_REGISTER_FP_TYPE(tensorwrapper::types::afloat);
 WTF_REGISTER_FP_TYPE(tensorwrapper::types::adouble);
 WTF_REGISTER_FP_TYPE(tensorwrapper::types::tafloat);
 WTF_REGISTER_FP_TYPE(tensorwrapper::types::tadouble);
+WTF_REGISTER_FP_TYPE(tensorwrapper::types::tmfloat);
+WTF_REGISTER_FP_TYPE(tensorwrapper::types::tmdouble);
 
 #else
 template<typename T>
@@ -194,6 +228,10 @@ template<typename T>
 using thresholded_affine_type = T;
 using tafloat                 = float;
 using tadouble                = double;
+template<typename T>
+using taylor_model_type = T;
+using tmfloat           = float;
+using tmdouble          = double;
 
 using floating_point_types = std::tuple<float, double>;
 
@@ -209,6 +247,9 @@ constexpr bool is_affine_v = false;
 template<typename T>
 constexpr bool is_thresholded_affine_v = false;
 
+template<typename T>
+constexpr bool is_taylor_model_v = false;
+
 template<typename T>
 constexpr bool is_uq_type_v = false;
 
@@ -227,6 +268,11 @@ T uq_upper(const T& value) {
     return value;
 }
 
+template<typename T>
+T uq_lower(const T& value) {
+    return value;
+}
+
 template<typename T, typename U>
 bool strictly_less(const T& lhs, const U& rhs) {
     return lhs < rhs;
diff --git a/cxx/src/tensorwrapper/buffer/detail_/hash_utilities.hpp b/cxx/src/tensorwrapper/buffer/detail_/hash_utilities.hpp
index 63d6b99a..236bc304 100644
--- a/cxx/src/tensorwrapper/buffer/detail_/hash_utilities.hpp
+++ b/cxx/src/tensorwrapper/buffer/detail_/hash_utilities.hpp
@@ -16,6 +16,7 @@
 
 #pragma once
 #include <boost/container_hash/hash.hpp>
+#include <string>
 #include <tensorwrapper/types/floating_point.hpp>
 
 /** @namespace tensorwrapper::buffer::detail_::hash_utilities
@@ -104,6 +105,43 @@ void hash_input(hash_type& seed,
     hash_input(seed, value.threshold());
 }
 
+/** @brief Specialization for sigma::TaylorModel values
+ *
+ *  @tparam T The floating point type of the model's coefficients
+ *  @param[in,out] seed The initial value of the hash, which is overwritten when
+ *                      the new value is added.
+ *  @param[in] value The new TaylorModel value being hashed and combined with
+ *                   the seed.
+ *
+ *  @return The updated hash value
+ *
+ *  @throw none No throw guarantee
+ */
+template<typename T>
+void hash_input(hash_type& seed, const types::taylor_model_type<T>& value) {
+    // Unlike the other UQ types, TaylorModel has an explicit empty state
+    // (representing the empty set, distinct from a zero-width [0,0]
+    // remainder); remainder().lower()/upper() throw std::domain_error for an
+    // empty model, which would violate this function's no-throw guarantee.
+    if(value.empty()) {
+        hash_input(seed, std::string("empty_taylor_model"));
+        return;
+    }
+    hash_input(seed, value.constant());
+    // coefficients() is a std::map keyed on Monomial (which orders variables
+    // lexicographically), so iteration order is already stable -- unlike the
+    // unordered_map cases above, no order-independent folding is needed.
+    for(const auto& [mono, coeff] : value.coefficients()) {
+        for(const auto& [var, exponent] : mono.exponents()) {
+            hash_input(seed, var);
+            hash_input(seed, exponent);
+        }
+        hash_input(seed, coeff);
+    }
+    hash_input(seed, value.remainder().lower());
+    hash_input(seed, value.remainder().upper());
+}
+
 #endif
 
 class HashVisitor {
diff --git a/tests/cxx/unit_tests/tensorwrapper/buffer/contiguous.cpp b/tests/cxx/unit_tests/tensorwrapper/buffer/contiguous.cpp
index b3276f7e..e4f0004c 100644
--- a/tests/cxx/unit_tests/tensorwrapper/buffer/contiguous.cpp
+++ b/tests/cxx/unit_tests/tensorwrapper/buffer/contiguous.cpp
@@ -907,3 +907,34 @@ TEST_CASE("interval contraction") {
     }
 #endif
 }
+
+TEST_CASE("taylor model contraction") {
+#ifdef ENABLE_SIGMA
+    using tm_type = sigma::TaylorModel<double>;
+
+    // Matrix-vector product with the identity, so the correct result is
+    // exactly a_buffer regardless of how TaylorModel's polynomial/remainder
+    // machinery represents each intermediate product.
+    tensorwrapper::shape::Smooth a_shape({2});
+    tensorwrapper::shape::Smooth b_shape({2, 2});
+
+    std::vector a_buffer{tm_type(1.0), tm_type(2.0)};
+    std::vector b_buffer{tm_type(1.0), tm_type(0.0), tm_type(0.0),
+                         tm_type(1.0)};
+
+    buffer::Contiguous a_tensor(a_buffer, a_shape);
+    buffer::Contiguous b_tensor(b_buffer, b_shape);
+
+    std::vector<tm_type> c_buffer(a_shape.size(), tm_type(0.0));
+    buffer::Contiguous c_tensor(c_buffer, a_shape);
+    c_tensor("k") = a_tensor("i") * b_tensor("k,i");
+
+    std::vector<double> corr{1.0, 2.0};
+
+    auto c_tensor_buffer = buffer::get_raw_data<tm_type>(c_tensor);
+    REQUIRE(c_tensor_buffer.size() == corr.size());
+    for(size_t i = 0; i < c_tensor_buffer.size(); ++i) {
+        REQUIRE(c_tensor_buffer[i].contains(corr[i]));
+    }
+#endif
+}
diff --git a/tests/cxx/unit_tests/tensorwrapper/buffer/detail_/hash_utilities.cpp b/tests/cxx/unit_tests/tensorwrapper/buffer/detail_/hash_utilities.cpp
index f7486b17..4fd4c334 100644
--- a/tests/cxx/unit_tests/tensorwrapper/buffer/detail_/hash_utilities.cpp
+++ b/tests/cxx/unit_tests/tensorwrapper/buffer/detail_/hash_utilities.cpp
@@ -75,6 +75,21 @@ TEMPLATE_LIST_TEST_CASE("hash_input", "", types::floating_point_types) {
         hash_input(corr, value.affine());
         hash_input(corr, value.threshold());
         REQUIRE(seed == corr);
+    } else if constexpr(types::is_taylor_model_v<TestType>) {
+        value_type value(1.0, 2.0);
+        hash_input(seed, value);
+        hash_type corr{0};
+        boost::hash_combine(corr, value.constant());
+        for(const auto& [mono, coeff] : value.coefficients()) {
+            for(const auto& [var, exponent] : mono.exponents()) {
+                boost::hash_combine(corr, var);
+                boost::hash_combine(corr, exponent);
+            }
+            boost::hash_combine(corr, coeff);
+        }
+        boost::hash_combine(corr, value.remainder().lower());
+        boost::hash_combine(corr, value.remainder().upper());
+        REQUIRE(seed == corr);
     } else if constexpr(types::is_uq_type_v<TestType>) {
         throw std::runtime_error("UQ type not registered for hash_input");
     } else {
diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp
index bc02e78f..df789480 100644
--- a/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp
+++ b/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp
@@ -71,7 +71,12 @@ TEMPLATE_LIST_TEST_CASE("add_noise", "", types::floating_point_types) {
         const double t = 0.01;
         auto out1      = add_noise<TestType>(matrix, t, 99);
         auto out2      = add_noise<TestType>(matrix, t, 99);
-        if constexpr(types::is_interval_v<TestType>) {
+        if constexpr(types::is_uq_type_v<TestType>) {
+            // Two independent calls create values with the same observable
+            // bounds but different internal error-source/error-symbol IDs,
+            // so compare bounds via the type-generic uq_* accessors (which
+            // together fully determine each UQ type's observable state)
+            // rather than the raw values or their internal representations.
             using wtf::fp::float_cast;
             auto b1 = make_contiguous(out1.buffer());
             auto b2 = make_contiguous(out2.buffer());
@@ -79,37 +84,12 @@ TEMPLATE_LIST_TEST_CASE("add_noise", "", types::floating_point_types) {
                 for(std::size_t j = 0; j < 2; ++j) {
                     const auto v1 = float_cast<TestType>(b1.get_elem({i, j}));
                     const auto v2 = float_cast<TestType>(b2.get_elem({i, j}));
-                    REQUIRE(v1.lower() == Catch::Approx(v2.lower()));
-                    REQUIRE(v1.upper() == Catch::Approx(v2.upper()));
-                }
-            }
-        } else if constexpr(types::is_uncertain_v<TestType>) {
-            // Two independent calls create values with the same mean/sd but
-            // different internal error-source IDs; compare observables
-            // directly.
-            using wtf::fp::float_cast;
-            auto b1 = make_contiguous(out1.buffer());
-            auto b2 = make_contiguous(out2.buffer());
-            for(std::size_t i = 0; i < 2; ++i) {
-                for(std::size_t j = 0; j < 2; ++j) {
-                    const auto v1 = float_cast<TestType>(b1.get_elem({i, j}));
-                    const auto v2 = float_cast<TestType>(b2.get_elem({i, j}));
-                    REQUIRE(v1.mean() == Catch::Approx(v2.mean()));
-                    REQUIRE(v1.sd() == Catch::Approx(v2.sd()));
-                }
-            }
-        } else if constexpr(types::is_affine_v<TestType> ||
-                            types::is_thresholded_affine_v<TestType>) {
-            // Same reason as uncertain: independent error-symbol IDs differ.
-            using wtf::fp::float_cast;
-            auto b1 = make_contiguous(out1.buffer());
-            auto b2 = make_contiguous(out2.buffer());
-            for(std::size_t i = 0; i < 2; ++i) {
-                for(std::size_t j = 0; j < 2; ++j) {
-                    const auto v1 = float_cast<TestType>(b1.get_elem({i, j}));
-                    const auto v2 = float_cast<TestType>(b2.get_elem({i, j}));
-                    REQUIRE(v1.center() == Catch::Approx(v2.center()));
-                    REQUIRE(v1.radius() == Catch::Approx(v2.radius()));
+                    REQUIRE(types::uq_center(v1) ==
+                            Catch::Approx(types::uq_center(v2)));
+                    REQUIRE(types::uq_lower(v1) ==
+                            Catch::Approx(types::uq_lower(v2)));
+                    REQUIRE(types::uq_upper(v1) ==
+                            Catch::Approx(types::uq_upper(v2)));
                 }
             }
         } else {
diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp
index 356d1792..3e041aab 100644
--- a/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp
+++ b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp
@@ -33,7 +33,7 @@ template<typename T>
 constexpr double eigen_system_tol =
   std::is_same_v<T, float> || std::is_same_v<T, types::ufloat> ||
       std::is_same_v<T, types::ifloat> || std::is_same_v<T, types::afloat> ||
-      std::is_same_v<T, types::tafloat> ?
+      std::is_same_v<T, types::tafloat> || std::is_same_v<T, types::tmfloat> ?
     1e-5 :
     1e-12;
 } // namespace
diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp
index 17cb1826..4bb2e46f 100644
--- a/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp
+++ b/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp
@@ -52,7 +52,8 @@ TEMPLATE_LIST_TEST_CASE("random_orthogonal_matrix", "",
                                    std::is_same_v<TestType, types::ufloat> ||
                                    std::is_same_v<TestType, types::ifloat> ||
                                    std::is_same_v<TestType, types::afloat> ||
-                                   std::is_same_v<TestType, types::tafloat> ?
+                                   std::is_same_v<TestType, types::tafloat> ||
+                                   std::is_same_v<TestType, types::tmfloat> ?
                                  1e-5 :
                                  1e-12;
         REQUIRE(approximately_equal(product, ident, tol));
diff --git a/tests/cxx/unit_tests/tensorwrapper/testing/helpers.hpp b/tests/cxx/unit_tests/tensorwrapper/testing/helpers.hpp
index d2db8ad2..4d4b7283 100644
--- a/tests/cxx/unit_tests/tensorwrapper/testing/helpers.hpp
+++ b/tests/cxx/unit_tests/tensorwrapper/testing/helpers.hpp
@@ -53,7 +53,8 @@ bool elements_equal(const T& lhs, const U& rhs) {
         }
     }();
     if constexpr(types::is_affine_v<concrete_t> ||
-                 types::is_thresholded_affine_v<concrete_t>) {
+                 types::is_thresholded_affine_v<concrete_t> ||
+                 types::is_taylor_model_v<concrete_t>) {
         return lv.range() == rhs.range();
     } else {
         return lv == rhs;
@@ -62,7 +63,7 @@ bool elements_equal(const T& lhs, const U& rhs) {
 
 template<typename T>
 constexpr double default_tolerance() {
-    if constexpr(types::is_affine_v<T> || types::is_thresholded_affine_v<T>) {
+    if constexpr(types::is_uq_type_v<T>) {
         // pow() is implemented as exp(log(x)*n); float-precision accumulates
         // ~1e-3 absolute error for values like 42^2 = 1764.
         return 1e-3;