diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 3fe48c1b021..ca69aae110a 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -644,8 +644,24 @@ void Tally::set_scores(const vector<std::string>& scores)
       break;
 
     case HEATING:
-      if (settings::photon_transport)
-        estimator_ = TallyEstimator::COLLISION;
+      if (settings::photon_transport) {
+        // Photon heating requires a collision estimator (analog energy
+        // balance). However, if the tally only scores neutrons, we can
+        // keep the tracklength estimator since neutron heating uses
+        // kerma coefficients which support tracklength scoring.
+        bool neutron_only = false;
+        for (auto i_filt : filters_) {
+          auto pf = dynamic_cast<ParticleFilter*>(
+            model::tally_filters[i_filt].get());
+          if (pf && pf->particles().size() == 1 &&
+              pf->particles()[0].is_neutron()) {
+            neutron_only = true;
+            break;
+          }
+        }
+        if (!neutron_only)
+          estimator_ = TallyEstimator::COLLISION;
+      }
       break;
 
     case SCORE_PULSE_HEIGHT: {
diff --git a/tests/unit_tests/test_tally.py b/tests/unit_tests/test_tally.py
index 52647a36c7c..ab6525794ef 100644
--- a/tests/unit_tests/test_tally.py
+++ b/tests/unit_tests/test_tally.py
@@ -17,3 +17,47 @@ def test_tally_init_args():
     assert tally.filters == [filter]
     assert tally.nuclides == ['U235']
     assert tally.estimator == 'tracklength'
+
+
+def test_heating_estimator_with_photon_transport(run_in_tmpdir):
+    """Test that a neutron-only heating tally uses the tracklength estimator
+    even when photon transport is enabled.
+
+    Without a neutron particle filter, photon heating requires the collision
+    estimator (analog energy balance). But neutron heating has kerma
+    coefficients that support tracklength scoring, so a neutron-only tally
+    should keep the tracklength estimator.
+    """
+    mat = openmc.Material()
+    mat.add_nuclide('Si28', 1.0)
+    mat.set_density('g/cm3', 2.3)
+    sph = openmc.Sphere(r=10.0, boundary_type='vacuum')
+    cell = openmc.Cell(fill=mat, region=-sph)
+
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell])
+    model.settings.run_mode = 'fixed source'
+    model.settings.particles = 100
+    model.settings.batches = 1
+    model.settings.photon_transport = True
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point()
+    )
+
+    neutron_filter = openmc.ParticleFilter('neutron')
+
+    # Neutron-only heating — should get tracklength estimator
+    t_neutron = openmc.Tally(name='heating_neutron')
+    t_neutron.filters = [neutron_filter]
+    t_neutron.scores = ['heating']
+
+    # No particle filter — should get collision estimator
+    t_all = openmc.Tally(name='heating_all')
+    t_all.scores = ['heating']
+
+    model.tallies = [t_neutron, t_all]
+
+    sp_filename = model.run()
+    with openmc.StatePoint(sp_filename) as sp:
+        assert sp.tallies[t_neutron.id].estimator == 'tracklength'
+        assert sp.tallies[t_all.id].estimator == 'collision'