diff --git a/doc/changelog.rst b/doc/changelog.rst
index a505f3c4..0b450ead 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -25,6 +25,8 @@ General
    - Switched :ref:`mjd_inverseFD` to use the CSR-format ``mjData.M`` representation instead of the legacy ``mjData.qM``
      for the mass matrix derivative. This changes the shape of the ``DmDq`` parameter from ``(nv x nM)`` to
      ``(nv x nC)``.
+   - Removed the legacy sparse ancestor-walk inertia matrix ``mjData.qM``. The joint-space inertia matrix is now stored
+     exclusively in the compressed sparse row (CSR) format ``mjData.M``.
 
 Version 3.10.0 (June 22, 2026)
 ------------------------------
diff --git a/doc/includes/references.h b/doc/includes/references.h
index 11afb93a..9887de34 100644
--- a/doc/includes/references.h
+++ b/doc/includes/references.h
@@ -223,8 +223,7 @@ typedef struct mjData_ {
 
   // computed by mj_fwdPosition/mj_makeM
   mjtNum* crb;               // com-based composite inertia and mass             (nbody x 10)
-  mjtNum* qM;                // inertia (sparse)                                 (nM x 1)
-  mjtNum* M;                 // reduced inertia (compressed sparse row)          (nC x 1)
+  mjtNum* M;                 // inertia (sparse)                                 (nC x 1)
 
   // computed by mj_fwdPosition/mj_factorM
   mjtNum* qLD;               // L'*D*L factorization of M (sparse)               (nC x 1)
@@ -273,7 +272,7 @@ typedef struct mjData_ {
   mjtNum* qDeriv;            // d (passive + actuator - bias) / d qvel           (nD x 1)
 
   // computed by mj_implicit/mju_factorLUSparse
-  mjtNum* qLU;               // sparse LU of (qM - dt*qDeriv)                    (nD x 1)
+  mjtNum* qLU;               // sparse LU of (M - dt*qDeriv)                     (nD x 1)
 
   //-------------------- POSITION, VELOCITY, CONTROL/ACCELERATION dependent
 
diff --git a/include/mujoco/mjdata.h b/include/mujoco/mjdata.h
index d95d4e5a..9c9e9938 100644
--- a/include/mujoco/mjdata.h
+++ b/include/mujoco/mjdata.h
@@ -247,8 +247,7 @@ typedef struct mjData_ {
 
   // computed by mj_fwdPosition/mj_makeM
   mjtNum* crb;               // com-based composite inertia and mass             (nbody x 10)
-  mjtNum* qM;                // inertia (sparse)                                 (nM x 1)
-  mjtNum* M;                 // reduced inertia (compressed sparse row)          (nC x 1)
+  mjtNum* M;                 // inertia (sparse)                                 (nC x 1)
 
   // computed by mj_fwdPosition/mj_factorM
   mjtNum* qLD;               // L'*D*L factorization of M (sparse)               (nC x 1)
@@ -297,7 +296,7 @@ typedef struct mjData_ {
   mjtNum* qDeriv;            // d (passive + actuator - bias) / d qvel           (nD x 1)
 
   // computed by mj_implicit/mju_factorLUSparse
-  mjtNum* qLU;               // sparse LU of (qM - dt*qDeriv)                    (nD x 1)
+  mjtNum* qLU;               // sparse LU of (M - dt*qDeriv)                     (nD x 1)
 
   //-------------------- POSITION, VELOCITY, CONTROL/ACCELERATION dependent
 
diff --git a/include/mujoco/mjxmacro.h b/include/mujoco/mjxmacro.h
index fbd34918..620dcf92 100644
--- a/include/mujoco/mjxmacro.h
+++ b/include/mujoco/mjxmacro.h
@@ -875,7 +875,6 @@
     X   ( int,       moment_colind,     nJmom,       1           ) \
     X   ( mjtNum,    actuator_moment,   nJmom,       1           ) \
     XNV ( mjtNum,    crb,               nbody,       10          ) \
-    XNV ( mjtNum,    qM,                nM,          1           ) \
     XNV ( mjtNum,    M,                 nC,          1           ) \
     XNV ( mjtNum,    qLD,               nC,          1           ) \
     X   ( mjtNum,    qLDiagInv,         nv,          1           ) \
diff --git a/mjx/mujoco/mjx/_src/io.py b/mjx/mujoco/mjx/_src/io.py
index 13df4f7f..3cc678dc 100644
--- a/mjx/mujoco/mjx/_src/io.py
+++ b/mjx/mujoco/mjx/_src/io.py
@@ -1297,12 +1297,11 @@ def _get_data_into_warp(
         value = value.reshape((-1, 9))
       # elif field.name == 'efc_J':  # TODO(btaba): add this back
       # elif field.name.startswith('efc_'):  # TODO(btaba): add this back
-      # TODO(btaba): qM, qLD, qLDiagInv
+      # TODO(btaba): qLD, qLDiagInv
 
       if field.name in (
           'actuator_moment',
           'contact',
-          'qM',
           'qLD',
           'qLU',
           'qLDiagInv',
@@ -1472,10 +1471,6 @@ def _get_data_into(
       else:
         setattr(result_i, field.name, value)
 
-    if hasattr(result_i, 'qM'):
-      result_i.qM.fill(0.0)
-      result_i.qM[m.mapM2M] = result_i.M
-
     # recalculate qLD and qLDiagInv as MJX and MuJoCo have different
     # representations of the Cholesky decomposition.
     mujoco.mj_factorM(m, result_i)
diff --git a/mjx/mujoco/mjx/_src/io_test.py b/mjx/mujoco/mjx/_src/io_test.py
index b7ca1257..7a21f911 100644
--- a/mjx/mujoco/mjx/_src/io_test.py
+++ b/mjx/mujoco/mjx/_src/io_test.py
@@ -671,7 +671,6 @@ class DataIOTest(parameterized.TestCase):
     np.testing.assert_allclose(d_2.xpos, d.xpos)
     np.testing.assert_allclose(d_2.cvel, d.cvel)
     np.testing.assert_allclose(d_2.cdof_dot, d.cdof_dot)
-    np.testing.assert_allclose(d_2.qM, d.qM)
     np.testing.assert_allclose(d_2.qLD, d.qLD, atol=1e-6)
     np.testing.assert_allclose(d_2.qLDiagInv, d.qLDiagInv, atol=1e-6)
 
@@ -782,7 +781,6 @@ class DataIOTest(parameterized.TestCase):
     # check a few fields
     np.testing.assert_allclose(d_2.qpos, d.qpos)
     np.testing.assert_allclose(d_2.xpos, d.xpos)
-    np.testing.assert_allclose(d_2.qM, d.qM)
 
     # only 1 contact active
     self.assertEqual(d_2.contact.dist.shape, (1,))
diff --git a/python/mujoco/introspect/structs.py b/python/mujoco/introspect/structs.py
index d00eb31f..f51c3926 100644
--- a/python/mujoco/introspect/structs.py
+++ b/python/mujoco/introspect/structs.py
@@ -6115,20 +6115,12 @@ STRUCTS: Mapping[str, StructDecl] = dict([
                  doc='com-based composite inertia and mass',
                  array_extent=('nbody', 10),
              ),
-             StructFieldDecl(
-                 name='qM',
-                 type=PointerType(
-                     inner_type=ValueType(name='mjtNum'),
-                 ),
-                 doc='inertia (sparse)',
-                 array_extent=('nM',),
-             ),
              StructFieldDecl(
                  name='M',
                  type=PointerType(
                      inner_type=ValueType(name='mjtNum'),
                  ),
-                 doc='reduced inertia (compressed sparse row)',
+                 doc='inertia (sparse)',
                  array_extent=('nC',),
              ),
              StructFieldDecl(
@@ -6328,7 +6320,7 @@ STRUCTS: Mapping[str, StructDecl] = dict([
                  type=PointerType(
                      inner_type=ValueType(name='mjtNum'),
                  ),
-                 doc='sparse LU of (qM - dt*qDeriv)',
+                 doc='sparse LU of (M - dt*qDeriv)',
                  array_extent=('nD',),
              ),
              StructFieldDecl(
diff --git a/src/engine/engine_core_smooth.c b/src/engine/engine_core_smooth.c
index 2d4ac47f..2f118b0d 100644
--- a/src/engine/engine_core_smooth.c
+++ b/src/engine/engine_core_smooth.c
@@ -1871,7 +1871,6 @@ void mj_makeM(const mjModel* m, mjData* d) {
   TM_START;
   mj_crb(m, d);
   mj_tendonArmature(m, d);
-  mju_scatter(d->qM, d->M, m->mapM2M, m->nC);  // TODO(tassa): scatter only awake dofs
   TM_END(mjTIMER_POS_INERTIA);
 }
 
diff --git a/src/engine/engine_inverse.c b/src/engine/engine_inverse.c
index 2bfbe3c0..bc10ea45 100644
--- a/src/engine/engine_inverse.c
+++ b/src/engine/engine_inverse.c
@@ -132,10 +132,10 @@ static void mj_discreteAcc(const mjModel* m, mjData* d) {
     // compute qDeriv
     mjd_smooth_vel(m, d, /* flg_bias = */ 1);
 
-    // gather qLU <- qM (lower to full)
+    // gather qLU <- M (lower to full)
     mju_gatherMasked(d->qLU, d->M, m->mapM2D, nD);
 
-    // set qLU = qM - dt*qDeriv
+    // set qLU = M - dt*qDeriv
     mju_addToScl(d->qLU, d->qDeriv, -m->opt.timestep, m->nD);
 
     // set qfrc = qLU * qacc
diff --git a/src/engine/engine_io.c b/src/engine/engine_io.c
index 0c455674..6f5139cf 100644
--- a/src/engine/engine_io.c
+++ b/src/engine/engine_io.c
@@ -1446,9 +1446,6 @@ static void _resetData(const mjModel* m, mjData* d, unsigned char debug_value) {
     }
   }
 
-  // zero out qM, special case because scattering from M skips simple body off-diagonals
-  mju_zero(d->qM, m->nM);
-
   // copy qpos0 from model
   if (m->qpos0) {
     mju_copy(d->qpos, m->qpos0, m->nq);
diff --git a/src/engine/engine_print.c b/src/engine/engine_print.c
index 3e5e04ed..df259f2a 100644
--- a/src/engine/engine_print.c
+++ b/src/engine/engine_print.c
@@ -242,41 +242,6 @@ static void printBlockArray(const char* str, const mjtNum* data, int nc,
 }
 
 
-// print sparse inertia-like matrix
-static void printInertia(const char* str, const mjtNum* mat, const mjModel* m,
-                         FILE* fp, const char* float_format) {
-  int nv = m->nv;
-  // if no data, or too many rows to be visually useful, return
-  if (!mat || !nv || nv > 300) {
-    return;
-  }
-
-  // get length of string produced by float_format
-  char test[100];
-  int len = snprintf(test, sizeof(test), float_format, 0.0);
-
-  fprintf(fp, "%s\n", str);
-
-  for (int i=0; i < nv; i++) {
-    fprintf(fp, " ");
-    int adr = (i == nv-1) ? m->nM - 1 : m->dof_Madr[i+1] - 1;
-    for (int k=0; k <= i; k++) {
-      int j = i;
-      while (j != k && j >= 0) {
-        j = m->dof_parentid[j];
-      }
-      if (j == k) {
-        fprintf(fp, " ");
-        fprintf(fp, float_format, mat[adr--]);
-      } else {
-        for (int d=0; d < len+1; d++) fprintf(fp, " ");
-      }
-    }
-    fprintf(fp, "\n");
-  }
-  fprintf(fp, "\n");
-}
-
 
 // print sparse matrix structure
 void mj_printSparsity(const char* str, int nr, int nc, const int* rowadr, const int* diag,
@@ -1464,7 +1429,6 @@ void mj_printFormattedData(const mjModel* m, const mjData* d, const char* filena
   printSparse("ACTUATOR_MOMENT", d->actuator_moment, m->nu, d->moment_rownnz,
               d->moment_rowadr, d->moment_colind, fp, float_format);
   printArray2d("CRB", m->nbody, 10, d->crb, fp, float_format);
-  printInertia("QM", d->qM, m, fp, float_format);
   printSparse("M", d->M, m->nv, m->M_rownnz,
               m->M_rowadr, m->M_colind, fp, float_format);
   printSparse("QLD", d->qLD, m->nv, m->M_rownnz,
diff --git a/test/engine/engine_core_smooth_test.cc b/test/engine/engine_core_smooth_test.cc
index 6cff4a3a..85e68955 100644
--- a/test/engine/engine_core_smooth_test.cc
+++ b/test/engine/engine_core_smooth_test.cc
@@ -227,7 +227,6 @@ TEST_F(CoreSmoothTest, TendonArmature) {
 
     // put only CRB inertia in M2
     mj_crb(m, d);
-    mju_scatter(d->qM, d->M, m->mapM2M, m->nC);
     vector<mjtNum> M2(nv * nv);
     mj_fullM(m, d, M2.data());
 
diff --git a/test/user/user_flex_test.cc b/test/user/user_flex_test.cc
index b1ae4459..cecbd4ee 100644
--- a/test/user/user_flex_test.cc
+++ b/test/user/user_flex_test.cc
@@ -390,9 +390,9 @@ TEST_F(UserFlexTest, TrilinearInterpolation) {
     EXPECT_NEAR(d1->flexvert_xpos[i], d2->flexvert_xpos[i], 1e-7);
   }
 
-  EXPECT_EQ(m1->nM, m2->nM);
-  for (int i = 0; i < m1->nM; ++i) {
-    EXPECT_EQ(d1->qM[i], d2->qM[i]);
+  EXPECT_EQ(m1->nC, m2->nC);
+  for (int i = 0; i < m1->nC; ++i) {
+    EXPECT_EQ(d1->M[i], d2->M[i]);
   }
 
   EXPECT_EQ(m1->nbody, m2->nbody);
diff --git a/unity/Runtime/Bindings/MjBindings.cs b/unity/Runtime/Bindings/MjBindings.cs
index e3f10d13..4adfa74d 100644
--- a/unity/Runtime/Bindings/MjBindings.cs
+++ b/unity/Runtime/Bindings/MjBindings.cs
@@ -5757,7 +5757,6 @@ public unsafe struct mjData_ {
   public int* moment_colind;
   public double* actuator_moment;
   public double* crb;
-  public double* qM;
   public double* M;
   public double* qLD;
   public double* qLDiagInv;
diff --git a/wasm/codegen/generated/bindings.cc b/wasm/codegen/generated/bindings.cc
index 88825c39..cb62406d 100644
--- a/wasm/codegen/generated/bindings.cc
+++ b/wasm/codegen/generated/bindings.cc
@@ -4664,7 +4664,6 @@ EMSCRIPTEN_BINDINGS(mujoco_bindings) {
     .property("qLD", &MjData::qLD)
     .property("qLDiagInv", &MjData::qLDiagInv)
     .property("qLU", &MjData::qLU)
-    .property("qM", &MjData::qM)
     .property("qacc", &MjData::qacc)
     .property("qacc_smooth", &MjData::qacc_smooth)
     .property("qacc_warmstart", &MjData::qacc_warmstart)
diff --git a/wasm/codegen/generated/bindings.h b/wasm/codegen/generated/bindings.h
index d0ad5a6f..05671f1d 100644
--- a/wasm/codegen/generated/bindings.h
+++ b/wasm/codegen/generated/bindings.h
@@ -6875,9 +6875,6 @@ struct MjData {
   emscripten::val crb() const {
     return emscripten::val(emscripten::typed_memory_view(model->nbody * 10, ptr_->crb));
   }
-  emscripten::val qM() const {
-    return emscripten::val(emscripten::typed_memory_view(model->nM, ptr_->qM));
-  }
   emscripten::val M() const {
     return emscripten::val(emscripten::typed_memory_view(model->nC, ptr_->M));
   }