diff --git a/CHANGELOG.md b/CHANGELOG.md
index da51812d13..a9f85a7acb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,8 @@
 
 ## Bug fixes
 
+- Heat of mixing can now be built with quadratic and quartic particle profiles; their radial concentration gradients are evaluated analytically instead of applying a finite-volume gradient to a reconstructed profile. ([#5902](https://github.com/pybamm-team/PyBaMM/pull/5902))
+
 - `"intercalation kinetics": "Marcus"` now gives a positive current for a positive overpotential, as the other kinetics do, and converts the reorganization energy from eV with F/RT, as Marcus-Hush-Chidsey does. Previously the oxidation and reduction terms were swapped and the energy was used in mixed units, so models with Marcus kinetics failed to solve. ([#5872](https://github.com/pybamm-team/PyBaMM/pull/5872))
 - Evaluating a processed variable at a 0-d array time, e.g. `solution["Voltage [V]"](np.array(10.0))`, now returns the same value as passing the float rather than raising an `IndexError`. ([#5762](https://github.com/pybamm-team/PyBaMM/pull/5762))
 - Fixed the particle concentration of the `"quartic profile"` particle model in models that resolve the electrode (e.g. DFN), which was missing a factor of the particle radius on the concentration gradient term and diverged to unphysical values. ([#5893](https://github.com/pybamm-team/PyBaMM/pull/5893))
diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
index fa5e0b590f..2ca5e6116d 100644
--- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
+++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
@@ -390,6 +390,25 @@ def _heat_of_mixing(self, variables):
             F = pybamm.constants.F.value
             pi = np.pi
 
+            def radial_gradient_squared(concentration):
+                # Polynomial profiles already define concentration as a function of
+                # radius; use its analytic derivative instead of a finite-volume one.
+                radial_variable = next(
+                    (
+                        symbol
+                        for symbol in concentration.pre_order()
+                        if isinstance(symbol, pybamm.SpatialVariable)
+                        and symbol.coord_sys == "spherical polar"
+                    ),
+                    None,
+                )
+                if radial_variable is not None:
+                    derivative = concentration.diff(radial_variable)
+                    return derivative**2
+
+                gradient = pybamm.grad(concentration)
+                return pybamm.inner(gradient, gradient)
+
             # Compute heat of mixing in negative electrode
             if self.options.electrode_types["negative"] == "planar":
                 Q_mix_s_n = pybamm.FullBroadcast(
@@ -408,7 +427,7 @@ def _heat_of_mixing(self, variables):
                     c_n = variables["Negative particle concentration [mol.m-3]"]
                     T_n = variables["Negative electrode temperature [K]"]
                 T_n_part = pybamm.PrimaryBroadcast(T_n, ["negative particle"])
-                dc_n_dr2 = pybamm.inner(pybamm.grad(c_n), pybamm.grad(c_n))
+                dc_n_dr2 = radial_gradient_squared(c_n)
                 D_n = self.param.n.prim.D(c_n, T_n_part)
                 dUeq_n = self.param.n.prim.U(
                     c_n / self.param.n.prim.c_max, T_n_part
@@ -431,7 +450,7 @@ def _heat_of_mixing(self, variables):
                 c_p = variables["Positive particle concentration [mol.m-3]"]
                 T_p = variables["Positive electrode temperature [K]"]
             T_p_part = pybamm.PrimaryBroadcast(T_p, ["positive particle"])
-            dc_p_dr2 = pybamm.inner(pybamm.grad(c_p), pybamm.grad(c_p))
+            dc_p_dr2 = radial_gradient_squared(c_p)
             D_p = self.param.p.prim.D(c_p, T_p_part)
             dUeq_p = self.param.p.prim.U(c_p / self.param.p.prim.c_max, T_p_part).diff(
                 c_p
diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
index 7f2609c4e8..ccfabd2aa7 100644
--- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
+++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
@@ -222,3 +222,25 @@ def R_contact(T):
             # cell heats up above the reference, so R (and hence the contact
             # overpotential magnitude) is larger than the constant case
             np.testing.assert_array_less(np.abs(dphi_const[1:]), np.abs(dphi_fn[1:]))
+
+    @pytest.mark.parametrize(
+        ("model_class", "particle_profile"),
+        [
+            (pybamm.lithium_ion.SPM, "Fickian diffusion"),
+            (pybamm.lithium_ion.SPM, "quadratic profile"),
+            (pybamm.lithium_ion.SPM, "quartic profile"),
+            (pybamm.lithium_ion.DFN, "Fickian diffusion"),
+            (pybamm.lithium_ion.DFN, "quadratic profile"),
+            (pybamm.lithium_ion.DFN, "quartic profile"),
+        ],
+    )
+    def test_heat_of_mixing_with_particle_profile(self, model_class, particle_profile):
+        model = model_class(
+            {
+                "particle": particle_profile,
+                "thermal": "lumped",
+                "heat of mixing": "true",
+            }
+        )
+
+        pybamm.Simulation(model).build()