diff --git a/README.md b/README.md
index b837f56..47f6b30 100644
--- a/README.md
+++ b/README.md
@@ -105,6 +105,8 @@ Static variants (`D1Scalar`–`D15Scalar`, `DD1Scalar`–`DD15Scalar`) avoid hea
 
 The convenience function `hj.variables(values, order=2)` automatically selects the appropriate static type when the number of variables is ≤ 15, and falls back to the dynamic variant otherwise.
 
+First-order types store only a value and a gradient. The Hessian accessors (`h`, `set_h`, `hm`, `set_hm`) therefore exist on second-order types only — in C++ they are constrained via `requires (order() == 2)`, and in Python they are absent from first-order classes.
+
 ### `SScalar` — Sparse dual numbers with named variables
 
 Stores first-order derivatives in a sparse map keyed by variable name (string). Useful when variables are identified by name rather than index, or when only a small subset of derivatives is non-zero.
diff --git a/include/hyperjet/hyperjet.h b/include/hyperjet/hyperjet.h
index 1ba906c..f748752 100644
--- a/include/hyperjet/hyperjet.h
+++ b/include/hyperjet/hyperjet.h
@@ -576,15 +576,23 @@ class DDScalar {
 
   void set_g(const index i, const Scalar value) { g(i) = value; }
 
-  auto &h(this auto &self, const index i) {
+  auto &h(this auto &self, const index i)
+    requires(order() == 2)
+  {
     assert(0 <= i && i < self.size() * (self.size() + 1) / 2);
 
     return self.m_data[1 + self.size() + i];
   }
 
-  void set_h(const index i, const Scalar value) { h(i) = value; }
+  void set_h(const index i, const Scalar value)
+    requires(order() == 2)
+  {
+    h(i) = value;
+  }
 
-  auto &h(this auto &self, const index i, const index j) {
+  auto &h(this auto &self, const index i, const index j)
+    requires(order() == 2)
+  {
     assert(0 <= i && i < self.size());
     assert(0 <= j && j < self.size());
 
@@ -597,7 +605,9 @@ class DDScalar {
     }
   }
 
-  void set_h(const index i, const index j, const Scalar value) {
+  void set_h(const index i, const index j, const Scalar value)
+    requires(order() == 2)
+  {
     h(i, j) = value;
   }
 
@@ -612,7 +622,9 @@ class DDScalar {
 
   Eigen::Ref<Vector> ag() { return Eigen::Map<Vector>(ptr() + 1, size()); }
 
-  Matrix hm(const std::string mode) const {
+  Matrix hm(const std::string mode) const
+    requires(order() == 2)
+  {
     Matrix result(size(), size());
 
     hm(mode, result);
@@ -620,7 +632,9 @@ class DDScalar {
     return result;
   }
 
-  void hm(const std::string mode, Eigen::Ref<Matrix> out) const {
+  void hm(const std::string mode, Eigen::Ref<Matrix> out) const
+    requires(order() == 2)
+  {
     index it = 0;
 
     for (index i = 0; i < size(); i++) {
@@ -646,7 +660,9 @@ class DDScalar {
     }
   }
 
-  void set_hm(const Eigen::Ref<const Matrix> &value) {
+  void set_hm(const Eigen::Ref<const Matrix> &value)
+    requires(order() == 2)
+  {
     index it = 0;
 
     for (index i = 0; i < size(); i++) {
diff --git a/python/src/common.h b/python/src/common.h
index 827c76c..789e17b 100644
--- a/python/src/common.h
+++ b/python/src/common.h
@@ -92,19 +92,23 @@ template <typename T> auto bind(py::module &m, const std::string &name) {
       .def("__pow__", &T::pow)
       .def("__repr__", &T::to_string)
       .def("abs", &T::abs)
-      .def("eval", &T::eval, "d"_a)
-      .def(
-          "h",
-          [](const T &self, hj::index row, hj::index col) ->
-          typename T::Scalar { return self.h(row, col); },
-          "row"_a, "col"_a)
-      .def("set_h",
-           py::overload_cast<hj::index, hj::index, typename T::Scalar>(
-               &T::set_h),
-           "row"_a, "col"_a, "value"_a)
-      .def("hm", py::overload_cast<std::string>(&T::hm, py::const_),
-           "mode"_a = "full")
-      .def("set_hm", &T::set_hm, "value"_a);
+      .def("eval", &T::eval, "d"_a);
+
+  // methods: Hessian (second order only)
+  if constexpr (T::order() == 2) {
+    cls.def(
+           "h",
+           [](const T &self, hj::index row, hj::index col) ->
+           typename T::Scalar { return self.h(row, col); },
+           "row"_a, "col"_a)
+        .def("set_h",
+             py::overload_cast<hj::index, hj::index, typename T::Scalar>(
+                 &T::set_h),
+             "row"_a, "col"_a, "value"_a)
+        .def("hm", py::overload_cast<std::string>(&T::hm, py::const_),
+             "mode"_a = "full")
+        .def("set_hm", &T::set_hm, "value"_a);
+  }
 
   if constexpr (T::is_dynamic()) {
     cls.def("resize", &T::resize, "size"_a)
diff --git a/python/tests/test_DDScalar.py b/python/tests/test_DDScalar.py
index e866bd1..1bf099d 100644
--- a/python/tests/test_DDScalar.py
+++ b/python/tests/test_DDScalar.py
@@ -543,6 +543,28 @@ def test_ndarray(ctx):
     ctx.check(u, [3, 4, 5, 6, 7, 8])
 
 
+@pytest.mark.parametrize("ctx", **test_data)
+def test_hessian_access(ctx):
+    u = ctx.from_data([1, 2, 3, 4, 5, 6])
+
+    if ctx.dtype.order == 2:
+        assert_equal(u.h(0, 0), 4)
+        assert_equal(u.hm(), [[4, 5], [5, 6]])
+        return
+
+    with pytest.raises(AttributeError):
+        u.h(0, 0)
+
+    with pytest.raises(AttributeError):
+        u.set_h(0, 0, 1)
+
+    with pytest.raises(AttributeError):
+        u.hm()
+
+    with pytest.raises(AttributeError):
+        u.set_hm([[1, 0], [0, 1]])
+
+
 def test_is_dynamic():
     assert_equal(static_set_2.u1.is_dynamic, False)
     assert_equal(dynamic_set_2.u1.is_dynamic, True)
@@ -1220,6 +1242,16 @@ def test_dd(ctx):
     assert_equal(dd, v.hm())
 
 
+@pytest.mark.parametrize("ctx", **test_data)
+def test_dd_of_first_order(ctx):
+    if ctx.dtype.order == 2:
+        return
+
+    u = [ctx.u1, ctx.u2]
+
+    assert_equal(hj.dd(u), np.empty((2, 0, 0)))
+
+
 def test_f_of_scalar():
     assert_equal(hj.f(1), 1)
 
diff --git a/test/src/test.cpp b/test/src/test.cpp
index b316f37..92888bf 100644
--- a/test/src/test.cpp
+++ b/test/src/test.cpp
@@ -473,6 +473,53 @@ TEST_CASE("Norm") {
   REQUIRE(r.h(2, 2) == doctest::Approx(9.2320222391647280));
 }
 
+// The Hessian accessors index into m_data behind the gradient. First-order
+// scalars have no such storage, so the accessors must not exist for them.
+
+template <typename T>
+concept HasHessianEntry = requires(T a) { a.h(hyperjet::index(0)); };
+
+template <typename T>
+concept HasHessianElement =
+    requires(T a) { a.h(hyperjet::index(0), hyperjet::index(0)); };
+
+template <typename T>
+concept HasSetHessianEntry =
+    requires(T a) { a.set_h(hyperjet::index(0), typename T::Scalar(0)); };
+
+template <typename T>
+concept HasSetHessianElement = requires(T a) {
+  a.set_h(hyperjet::index(0), hyperjet::index(0), typename T::Scalar(0));
+};
+
+template <typename T>
+concept HasHessianMatrix = requires(const T a) { a.hm(std::string("full")); };
+
+template <typename T>
+concept HasSetHessianMatrix = requires(T a) { a.set_hm(typename T::Matrix()); };
+
+TEST_CASE("Hessian access is restricted to second order") {
+  CHECK(HasHessianEntry<DDScalar<2, double, 3>>);
+  CHECK(HasHessianElement<DDScalar<2, double, 3>>);
+  CHECK(HasSetHessianEntry<DDScalar<2, double, 3>>);
+  CHECK(HasSetHessianElement<DDScalar<2, double, 3>>);
+  CHECK(HasHessianMatrix<DDScalar<2, double, 3>>);
+  CHECK(HasSetHessianMatrix<DDScalar<2, double, 3>>);
+
+  CHECK(HasHessianElement<DDScalar<2, double, Dynamic>>);
+  CHECK(HasHessianMatrix<DDScalar<2, double, Dynamic>>);
+
+  CHECK_FALSE(HasHessianEntry<DDScalar<1, double, 3>>);
+  CHECK_FALSE(HasHessianElement<DDScalar<1, double, 3>>);
+  CHECK_FALSE(HasSetHessianEntry<DDScalar<1, double, 3>>);
+  CHECK_FALSE(HasSetHessianElement<DDScalar<1, double, 3>>);
+  CHECK_FALSE(HasHessianMatrix<DDScalar<1, double, 3>>);
+  CHECK_FALSE(HasSetHessianMatrix<DDScalar<1, double, 3>>);
+
+  CHECK_FALSE(HasHessianElement<DDScalar<1, double, Dynamic>>);
+  CHECK_FALSE(HasHessianMatrix<DDScalar<1, double, Dynamic>>);
+}
+
 const SScalar<double> s1(3.0, {{"x", 1.0}, {"y", 6.0}, {"z", 4.0}});
 const SScalar<double> s2(4.0, {{"x", 7.0}, {"y", 1.0}});
 const SScalar<double> s3(0.3, {{"x", 0.1}, {"y", 0.8}, {"z", 0.2}});