diff --git a/docs/Features/condensation_strategy_system.md b/docs/Features/condensation_strategy_system.md
index c2e24b88e..347556acb 100644
--- a/docs/Features/condensation_strategy_system.md
+++ b/docs/Features/condensation_strategy_system.md
@@ -12,12 +12,20 @@ into runnable pipelines. You can choose simultaneous isothermal updates or
 staggered two-pass Gauss-Seidel sweeps with theta-controlled first-pass
 fractions, batch partitioning, and gas-field updates for stability.
 
+CondensationLatentHeat mirrors the isothermal workflow for particle-resolved
+runs while adding latent-heat-aware mass transfer and per-step energy
+diagnostics. The step now supports a `dynamic_viscosity` override and records
+`last_latent_heat_energy` for parity with the isothermal API plus energy
+tracking.
+
 This feature is built around user-facing APIs exposed via `particula.dynamics`:
 
 - `CondensationStrategy` – abstract base defining `rate` / `step`.
 - `CondensationIsothermal` – simultaneous isothermal mass transfer.
 - `CondensationIsothermalStaggered` – two-pass staggered (Gauss-Seidel) update
   with theta modes (`"half"`, `"random"`, `"batch"`) and batching.
+- `CondensationLatentHeat` – latent-heat-corrected rate with per-step energy
+  diagnostics.
 - `CondensationIsothermalBuilder`, `CondensationIsothermalStaggeredBuilder` –
   fluent builders with validation and unit handling.
 - `CondensationFactory` – factory selecting a condensation strategy by name.
@@ -38,6 +46,8 @@ This feature is built around user-facing APIs exposed via `particula.dynamics`:
   fields are depleted.
 - **Pipeline-ready**: Use `MassCondensation` with `sub_steps` for tight coupling
   to other runnables in a single pipeline.
+- **Latent heat diagnostics**: Track per-step energy release for
+  particle-resolved runs with `CondensationLatentHeat`.
 
 ## Who It's For
 
@@ -65,6 +75,7 @@ par.dynamics.CondensationStrategy
 # Concrete implementations
 par.dynamics.CondensationIsothermal
 par.dynamics.CondensationIsothermalStaggered
+par.dynamics.CondensationLatentHeat
 
 # Builders and factory
 par.dynamics.CondensationIsothermalBuilder
@@ -131,6 +142,39 @@ particle, gas = iso.step(
 )
 ```
 
+### CondensationLatentHeat (energy diagnostics)
+
+`CondensationLatentHeat` mirrors the isothermal step but applies a latent-heat
+correction when a latent heat strategy (or scalar fallback) is provided. It
+records `last_latent_heat_energy` each step (positive for condensation,
+negative for evaporation) and accepts a `dynamic_viscosity` override for
+particle-resolved workflows.
+
+```python
+latent = par.dynamics.CondensationLatentHeat(
+    molar_mass=0.018,
+    diffusion_coefficient=2e-5,
+    accommodation_coefficient=1.0,
+    latent_heat=2.4e6,  # J/kg fallback
+)
+particle, gas = latent.step(
+    particle=particle,
+    gas_species=gas,
+    temperature=298.15,
+    pressure=101325.0,
+    time_step=1.0,
+    dynamic_viscosity=1.8e-5,
+)
+energy_released = latent.last_latent_heat_energy  # total energy [J]
+```
+
+When no latent heat strategy is configured (or a nonpositive scalar is
+provided), the step follows the isothermal path and reports
+`last_latent_heat_energy = 0.0`.
+
+`last_latent_heat_energy` records the total latent heat released per step
+(sum of dm × L), not an energy density.
+
 ### CondensationIsothermalStaggered (two-pass Gauss-Seidel)
 
 `CondensationIsothermalStaggered` splits each timestep into two passes. Theta
@@ -476,6 +520,8 @@ aerosol, gas = selective.step(
 | `accommodation_coefficient` | Mass accommodation coefficient (unitless). | `1.0` |
 | `update_gases` | Whether to deplete gas concentrations during step. | `True` |
 | `skip_partitioning_indices` | Species indices to exclude from partitioning. | `None` |
+| `latent_heat_strategy` | Optional latent heat strategy. | `None` |
+| `latent_heat` | Scalar fallback latent heat [J/kg]. | `0.0` |
 | `theta_mode` | Staggered theta selection: `"half"`, `"random"`, `"batch"`. | `"half"` (staggered) |
 | `num_batches` | Gauss-Seidel batch count (clipped to particle count). | `1` |
 | `shuffle_each_step` | Shuffle particle order each step (staggered). | `True` |
@@ -500,7 +546,7 @@ aerosol, gas = selective.step(
 
 - Staggered solver is Gauss-Seidel only; other solvers are not exposed.
 - Factory supports `"isothermal"` and `"isothermal_staggered"` only.
-- No latent-heat or temperature feedback; condensation is isothermal.
+- No temperature feedback; latent heat is diagnostic only.
 - Minimum-radius clamp (1e-10 m) enforces continuum validity; sub-continuum
   physics is out of scope.
 
diff --git a/docs/index.md b/docs/index.md
index abee15f70..2fea125dd 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -31,8 +31,8 @@ Whether you’re a researcher, educator, or industry expert, Particula is design
 - **Supporting non-isothermal condensation** with thermal resistance,
   latent-heat mass transfer rate utilities, latent-heat energy release
   bookkeeping, and the `CondensationLatentHeat` strategy with
-  latent-heat-corrected `mass_transfer_rate()`/`rate()` (step integration
-  pending).
+  latent-heat-corrected `mass_transfer_rate()`/`rate()` and an implemented
+  `step()` that tracks per-step latent heat energy.
 - **Interrogating your experimental data** to validate and expand your impact.
 - **Fostering open-source collaboration** to share ideas and build on each other’s work.
 
diff --git a/particula/dynamics/condensation/condensation_strategies.py b/particula/dynamics/condensation/condensation_strategies.py
index 387f53164..243a9ccb4 100644
--- a/particula/dynamics/condensation/condensation_strategies.py
+++ b/particula/dynamics/condensation/condensation_strategies.py
@@ -33,6 +33,7 @@
 
 from particula.dynamics.condensation.mass_transfer import (
     get_first_order_mass_transport_k,
+    get_latent_heat_energy_released,
     get_mass_transfer,
     get_mass_transfer_rate,
     get_mass_transfer_rate_latent_heat,
@@ -1761,7 +1762,9 @@ class CondensationLatentHeat(CondensationStrategy):
     Attributes:
         latent_heat_strategy_input: Strategy input provided at initialization.
         latent_heat_input: Raw latent heat value provided at initialization.
-        last_latent_heat_energy: Diagnostic latent heat energy tracker.
+        last_latent_heat_energy: Total latent heat energy released in the most
+            recent step [J]. Positive values indicate condensation and
+            negative values indicate evaporation. Overwritten each call.
     """
 
     # pylint: disable=R0913, R0917
@@ -1952,11 +1955,9 @@ def mass_transfer_rate(
         pressure_delta = self.calculate_pressure_delta(
             particle, gas_species, temperature, radius_with_fill
         )
-        if not np.all(np.isfinite(pressure_delta)):
-            raise ValueError(
-                "Non-finite pressure_delta computed for latent-heat "
-                "condensation."
-            )
+        pressure_delta = np.nan_to_num(
+            pressure_delta, posinf=0.0, neginf=0.0, nan=0.0
+        )
 
         if self._latent_heat_strategy is None:
             return get_mass_transfer_rate(
@@ -2036,6 +2037,144 @@ def rate(
         rates = self._apply_skip_partitioning(rates)
         return rates
 
+    def _update_latent_heat_energy(
+        self, mass_transfer: NDArray[np.float64], temperature: float
+    ) -> None:
+        """Store latent heat energy released for the current step.
+
+        Args:
+            mass_transfer: Per-particle mass transfer in kg for this step.
+            temperature: System temperature in Kelvin.
+        """
+        if self._latent_heat_strategy is None:
+            self.last_latent_heat_energy = 0.0
+            return
+
+        latent_heat = self._latent_heat_strategy.latent_heat(temperature)
+        self.last_latent_heat_energy = float(
+            np.sum(
+                get_latent_heat_energy_released(
+                    mass_transfer,
+                    latent_heat,
+                )
+            )
+        )
+
+    def _calculate_norm_conc(
+        self, volume: float, concentration: NDArray[np.float64]
+    ) -> NDArray[np.float64]:
+        """Validate inputs and compute normalized concentration.
+
+        Args:
+            volume: Particle volume for the single box [m^3].
+            concentration: Particle concentration array.
+
+        Returns:
+            Normalized concentration (concentration / volume).
+
+        Raises:
+            ValueError: If volume or concentration inputs are invalid.
+        """
+        if not np.isfinite(volume) or volume <= 0.0:
+            raise ValueError("volume must be finite and positive.")
+        if not np.all(np.isfinite(concentration)) or np.any(
+            concentration < 0.0
+        ):
+            raise ValueError(
+                "concentration must be finite and nonnegative for norm_conc."
+            )
+        return concentration / volume
+
+    def _normalize_mass_transfer_shape(
+        self,
+        mass_transfer: NDArray[np.float64],
+        species_mass: NDArray[np.float64],
+    ) -> NDArray[np.float64]:
+        """Normalize mass transfer shape to match species count.
+
+        Args:
+            mass_transfer: Per-particle mass transfer array.
+            species_mass: Particle mass array used to infer species count.
+
+        Returns:
+            Mass transfer array with consistent species dimension.
+        """
+        species_count = 1 if species_mass.ndim == 1 else species_mass.shape[1]
+        if mass_transfer.ndim == 1:
+            return mass_transfer.reshape(-1, species_count)
+        if mass_transfer.shape[1] > species_count:
+            return mass_transfer[:, :species_count]
+        if mass_transfer.shape[1] < species_count:
+            return np.broadcast_to(
+                mass_transfer, (mass_transfer.shape[0], species_count)
+            )
+        return mass_transfer
+
+    def _apply_mass_transfer_to_particles(
+        self,
+        particle_data: ParticleData,
+        mass_transfer: NDArray[np.float64],
+        norm_conc: NDArray[np.float64],
+    ) -> None:
+        """Apply mass transfer to particle data in place."""
+        if mass_transfer.ndim == 2:
+            nonzero_mask = norm_conc != 0
+            denom = norm_conc[:, np.newaxis]
+            where_mask = nonzero_mask[:, np.newaxis]
+        else:
+            denom = norm_conc
+            where_mask = norm_conc != 0
+        per_particle = np.divide(
+            mass_transfer,
+            denom,
+            out=np.zeros_like(mass_transfer),
+            where=where_mask,
+        )
+        np.add(
+            particle_data.masses[0],
+            per_particle,
+            out=particle_data.masses[0],
+        )
+        np.maximum(particle_data.masses[0], 0.0, out=particle_data.masses[0])
+
+    def _apply_mass_transfer_to_gas(
+        self, gas_data: GasData, mass_transfer: NDArray[np.float64]
+    ) -> None:
+        """Apply mass transfer to gas concentrations in place."""
+        np.subtract(
+            gas_data.concentration[0],
+            mass_transfer.sum(axis=0),
+            out=gas_data.concentration[0],
+        )
+        np.maximum(
+            gas_data.concentration[0],
+            0.0,
+            out=gas_data.concentration[0],
+        )
+
+    # pylint: disable=too-many-positional-arguments, too-many-arguments
+    @overload  # type: ignore[override]
+    def step(
+        self,
+        particle: ParticleRepresentation,
+        gas_species: GasSpecies,
+        temperature: float,
+        pressure: float,
+        time_step: float,
+        dynamic_viscosity: Optional[float] = None,
+    ) -> Tuple[ParticleRepresentation, GasSpecies]: ...
+
+    @overload
+    def step(
+        self,
+        particle: ParticleData,
+        gas_species: GasData,
+        temperature: float,
+        pressure: float,
+        time_step: float,
+        dynamic_viscosity: Optional[float] = None,
+    ) -> Tuple[ParticleData, GasData]: ...
+
     def step(
         self,
         particle: ParticleRepresentation | ParticleData,
@@ -2047,7 +2186,13 @@ def step(
     ) -> (
         Tuple[ParticleRepresentation, GasSpecies] | Tuple[ParticleData, GasData]
     ):
-        """Advance one condensation step (stub).
+        """Advance one condensation step with latent-heat diagnostics.
+
+        The mass transfer rate is computed, optional skip-partitioning applied,
+        and both the particle and gas states are updated while respecting
+        inventory limits. The per-step latent heat release is stored in
+        ``last_latent_heat_energy`` [J] (positive for condensation, negative
+        for evaporation) and overwritten on every call.
 
         Args:
             particle: Particle representation to update.
@@ -2055,15 +2200,69 @@ def step(
             temperature: System temperature in Kelvin.
             pressure: System pressure in Pascals.
             time_step: Integration timestep in seconds.
-            dynamic_viscosity: Optional dynamic viscosity override.
+            dynamic_viscosity: Optional dynamic viscosity override used in the
+                mass-transfer calculation.
 
         Returns:
             Tuple containing updated particle and gas species objects.
-
-        Raises:
-            NotImplementedError: Step is not implemented for latent-heat
-                condensation.
         """
-        raise NotImplementedError(
-            "CondensationLatentHeat.step is not implemented."
+        particle_data, particle_is_legacy = _unwrap_particle(particle)
+        gas_data, gas_is_legacy = _unwrap_gas(gas_species)
+        _require_matching_types(particle_is_legacy, gas_is_legacy)
+        _require_single_box(particle_data.n_boxes, "ParticleData")
+        _require_single_box(gas_data.n_boxes, "GasData")
+
+        mass_rate = self.mass_transfer_rate(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=temperature,
+            pressure=pressure,
+            dynamic_viscosity=dynamic_viscosity,
         )
+
+        mass_rate_array = np.atleast_1d(np.asarray(mass_rate, dtype=np.float64))
+
+        mass_rate_array = self._apply_skip_partitioning(mass_rate_array)
+
+        gas_mass_array: NDArray[np.float64] = np.atleast_1d(
+            np.asarray(gas_data.concentration[0], dtype=np.float64)
+        )
+        volume = particle_data.volume[0]
+        concentration = particle_data.concentration[0]
+        norm_conc = self._calculate_norm_conc(volume, concentration)
+        mass_transfer = get_mass_transfer(
+            mass_rate=mass_rate_array,
+            time_step=time_step,
+            gas_mass=gas_mass_array,
+            particle_mass=particle_data.masses[0],
+            particle_concentration=norm_conc,
+        )
+
+        mass_transfer = self._normalize_mass_transfer_shape(
+            mass_transfer,
+            particle_data.masses[0],
+        )
+
+        self._update_latent_heat_energy(mass_transfer, temperature)
+
+        if particle_is_legacy:
+            particle.add_mass(added_mass=mass_transfer)  # type: ignore[union-attr]
+        else:
+            self._apply_mass_transfer_to_particles(
+                particle_data,
+                mass_transfer,
+                norm_conc,
+            )
+        if self.update_gases:
+            if gas_is_legacy:
+                gas_species.add_concentration(  # type: ignore[union-attr]
+                    added_concentration=-mass_transfer.sum(axis=0)
+                )
+            else:
+                self._apply_mass_transfer_to_gas(gas_data, mass_transfer)
+        if particle_is_legacy:
+            return cast(
+                Tuple[ParticleRepresentation, GasSpecies],
+                (particle, gas_species),
+            )
+        return cast(Tuple[ParticleData, GasData], (particle_data, gas_data))
diff --git a/particula/dynamics/condensation/tests/condensation_strategies_test.py b/particula/dynamics/condensation/tests/condensation_strategies_test.py
index 0384c5764..dff15bd2f 100644
--- a/particula/dynamics/condensation/tests/condensation_strategies_test.py
+++ b/particula/dynamics/condensation/tests/condensation_strategies_test.py
@@ -22,6 +22,7 @@
     _unwrap_gas,
     _unwrap_particle,
 )
+from particula.dynamics.condensation.mass_transfer import get_mass_transfer
 from particula.gas.gas_data import GasData, from_species
 from particula.gas.latent_heat_strategies import ConstantLatentHeat
 from particula.gas.species import GasSpecies
@@ -1914,13 +1915,6 @@ def test_base_class_params_forwarded(self):
         self.assertEqual(cond.accommodation_coefficient, 0.7)
         self.assertFalse(cond.update_gases)
 
-    def test_stubs_raise_not_implemented(self):
-        """step() remains stubbed for latent heat strategy."""
-        cond = CondensationLatentHeat(molar_mass=0.018)
-
-        with self.assertRaises(NotImplementedError):
-            cond.step(None, None, 298.15, 101325, time_step=1.0)
-
     def test_mass_transfer_rate_isothermal_parity(self):
         """Latent heat strategy falls back to isothermal rates."""
         latent_cond = CondensationLatentHeat(
@@ -2132,7 +2126,7 @@ def test_rate_skip_partitioning_out_of_bounds_raises(self):
             )
 
     def test_nonfinite_pressure_delta_sanitized(self):
-        """Non-finite pressure deltas raise an error."""
+        """Non-finite pressure deltas are sanitized to zero."""
         gas_species, particle, _, _ = self._build_single_species_fixture(
             n_particles=3
         )
@@ -2143,13 +2137,50 @@ def test_nonfinite_pressure_delta_sanitized(self):
             "calculate_pressure_delta",
             return_value=np.array([np.nan, np.inf, -np.inf]),
         ):
-            with self.assertRaises(ValueError):
-                cond.mass_transfer_rate(
-                    particle=particle,
-                    gas_species=gas_species,
-                    temperature=self.temperature,
-                    pressure=self.pressure,
-                )
+            rates = cond.mass_transfer_rate(
+                particle=particle,
+                gas_species=gas_species,
+                temperature=self.temperature,
+                pressure=self.pressure,
+            )
+        self.assertTrue(np.all(np.isfinite(rates)))
+        self.assertTrue(np.allclose(rates, 0.0))
+
+    def test_nonfinite_pressure_delta_matches_isothermal(self):
+        """Latent heat sanitization mirrors isothermal behavior."""
+        gas_species, particle, _, _ = self._build_single_species_fixture(
+            n_particles=3
+        )
+        latent_cond = CondensationLatentHeat(molar_mass=self.molar_mass)
+        iso_cond = CondensationIsothermal(molar_mass=self.molar_mass)
+        nonfinite = np.array([np.nan, np.inf, -np.inf])
+
+        with (
+            patch.object(
+                CondensationLatentHeat,
+                "calculate_pressure_delta",
+                return_value=nonfinite,
+            ),
+            patch.object(
+                CondensationIsothermal,
+                "calculate_pressure_delta",
+                return_value=nonfinite,
+            ),
+        ):
+            latent_rates = latent_cond.mass_transfer_rate(
+                particle=particle,
+                gas_species=gas_species,
+                temperature=self.temperature,
+                pressure=self.pressure,
+            )
+            iso_rates = iso_cond.mass_transfer_rate(
+                particle=particle,
+                gas_species=gas_species,
+                temperature=self.temperature,
+                pressure=self.pressure,
+            )
+
+        np.testing.assert_allclose(latent_rates, iso_rates, rtol=1e-14)
 
     def test_nonfinite_latent_heat_raises(self):
         """Non-finite latent heat raises an error."""
@@ -2229,3 +2260,627 @@ def test_inherited_methods_work(self):
 
         self.assertGreater(mean_free_path, 0.0)
         self.assertGreater(knudsen_number, 0.0)
+
+    def test_step_mass_conservation_single_species(self):
+        """Single-species step conserves total mass."""
+        gas_species, particle, _, _ = self._build_single_species_fixture()
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+        )
+
+        initial_mass = (
+            particle.get_mass_concentration()
+            + gas_species.get_concentration().sum()
+        )
+        particle_new, gas_new = cond.step(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+        final_mass = (
+            particle_new.get_mass_concentration()
+            + gas_new.get_concentration().sum()
+        )
+
+        relative_error = abs(final_mass - initial_mass) / initial_mass
+        self.assertLess(relative_error, 1e-12)
+
+    def test_step_mass_conservation_multi_species(self):
+        """Multi-species step conserves total mass."""
+        particle = copy.deepcopy(self.particle)
+        gas_species = copy.deepcopy(self.gas_species)
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+        )
+
+        initial_mass = (
+            particle.get_mass_concentration()
+            + gas_species.get_concentration().sum()
+        )
+        particle_new, gas_new = cond.step(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+        final_mass = (
+            particle_new.get_mass_concentration()
+            + gas_new.get_concentration().sum()
+        )
+
+        relative_error = abs(final_mass - initial_mass) / initial_mass
+        self.assertLess(relative_error, 1e-12)
+
+    def test_step_mass_conservation_data_only(self):
+        """Data-only step conserves total mass."""
+        particle_data, gas_data = self._make_data_inputs()
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            activity_strategy=self.activity_strategy,
+            surface_strategy=self.surface_strategy,
+            vapor_pressure_strategy=self.vapor_pressure_strategy,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+        )
+
+        initial_mass = np.sum(
+            particle_data.masses[0]
+            * particle_data.concentration[0][:, np.newaxis]
+            / particle_data.volume[0]
+        ) + np.sum(gas_data.concentration[0])
+        particle_new, gas_new = cond.step(
+            particle=particle_data,
+            gas_species=gas_data,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+        final_mass = np.sum(
+            particle_new.masses[0]
+            * particle_new.concentration[0][:, np.newaxis]
+            / particle_new.volume[0]
+        ) + np.sum(gas_new.concentration[0])
+
+        relative_error = abs(final_mass - initial_mass) / initial_mass
+        self.assertLess(relative_error, 1e-12)
+
+    def test_step_rejects_nonpositive_volume(self):
+        """Step rejects nonpositive volumes for norm_conc."""
+        particle_data, gas_data = self._make_data_inputs()
+        particle_data.volume[0] = 0.0
+
+        with self.assertRaisesRegex(ValueError, "volume must be finite"):
+            self.data_strategy.step(
+                particle=particle_data,
+                gas_species=gas_data,
+                temperature=self.temperature,
+                pressure=self.pressure,
+                time_step=1.0,
+            )
+
+    def test_step_rejects_invalid_concentration_for_norm_conc(self):
+        """Step rejects non-finite or negative concentrations for norm_conc."""
+        particle_data, gas_data = self._make_data_inputs()
+
+        for bad_value in (-1.0, np.nan):
+            with self.subTest(bad_value=bad_value):
+                particle_copy = copy.deepcopy(particle_data)
+                particle_copy.concentration[0][0] = bad_value
+
+                with self.assertRaisesRegex(
+                    ValueError, "concentration must be finite"
+                ):
+                    self.data_strategy.step(
+                        particle=particle_copy,
+                        gas_species=gas_data,
+                        temperature=self.temperature,
+                        pressure=self.pressure,
+                        time_step=1.0,
+                    )
+
+    def test_step_energy_tracking_matches_dm_times_l(self):
+        """Latent heat energy matches dm × L for a step."""
+        particle = copy.deepcopy(self.particle)
+        gas_species = copy.deepcopy(self.gas_species)
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+        )
+
+        particle_data, _ = _unwrap_particle(particle)
+        gas_data, _ = _unwrap_gas(gas_species)
+        mass_rate = cond.mass_transfer_rate(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+        )
+        if isinstance(mass_rate, (int, float)):
+            mass_rate_array = np.array([mass_rate])
+        else:
+            mass_rate_array = mass_rate
+        mass_rate_array = cond._apply_skip_partitioning(mass_rate_array)
+        norm_conc = particle_data.concentration[0] / particle_data.volume[0]
+        mass_transfer = get_mass_transfer(
+            mass_rate=mass_rate_array,
+            time_step=1.0,
+            gas_mass=np.asarray(gas_data.concentration[0], dtype=np.float64),
+            particle_mass=particle_data.masses[0],
+            particle_concentration=norm_conc,
+        )
+
+        cond.step(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+        latent_heat = cond._latent_heat_strategy.latent_heat(self.temperature)
+
+        np.testing.assert_allclose(
+            cond.last_latent_heat_energy,
+            np.sum(mass_transfer * latent_heat),
+            rtol=1e-14,
+        )
+
+    def test_step_energy_tracking_update_gases_false(self):
+        """Energy tracking works when gas updates are disabled."""
+        gas_species, particle, _, _ = self._build_single_species_fixture()
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+            update_gases=False,
+        )
+
+        particle_data, _ = _unwrap_particle(particle)
+        gas_data, _ = _unwrap_gas(gas_species)
+        mass_rate = cond.mass_transfer_rate(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+        )
+        if isinstance(mass_rate, (int, float)):
+            mass_rate_array = np.array([mass_rate])
+        else:
+            mass_rate_array = mass_rate
+        mass_rate_array = cond._apply_skip_partitioning(mass_rate_array)
+        norm_conc = particle_data.concentration[0] / particle_data.volume[0]
+        mass_transfer = get_mass_transfer(
+            mass_rate=mass_rate_array,
+            time_step=1.0,
+            gas_mass=np.asarray(gas_data.concentration[0], dtype=np.float64),
+            particle_mass=particle_data.masses[0],
+            particle_concentration=norm_conc,
+        )
+        species_mass = particle_data.masses[0]
+        species_count = 1 if species_mass.ndim == 1 else species_mass.shape[1]
+        if mass_transfer.ndim == 1:
+            mass_transfer = mass_transfer.reshape(-1, species_count)
+        elif mass_transfer.shape[1] > species_count:
+            mass_transfer = mass_transfer[:, :species_count]
+        elif mass_transfer.shape[1] < species_count:
+            mass_transfer = np.broadcast_to(
+                mass_transfer, (mass_transfer.shape[0], species_count)
+            )
+        initial_conc = gas_species.get_concentration().copy()
+
+        cond.step(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+
+        np.testing.assert_allclose(
+            gas_species.get_concentration(), initial_conc, rtol=1e-15
+        )
+        latent_heat = cond._latent_heat_strategy.latent_heat(self.temperature)
+        np.testing.assert_allclose(
+            cond.last_latent_heat_energy,
+            np.sum(mass_transfer * latent_heat),
+            rtol=1e-14,
+        )
+
+    def test_step_energy_tracking_data_inputs(self):
+        """Latent heat energy tracks on data container inputs."""
+        particle_data, gas_data = self._make_data_inputs()
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            activity_strategy=self.activity_strategy,
+            surface_strategy=self.surface_strategy,
+            vapor_pressure_strategy=self.vapor_pressure_strategy,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+        )
+
+        mass_rate = cond.mass_transfer_rate(
+            particle=particle_data,
+            gas_species=gas_data,
+            temperature=self.temperature,
+            pressure=self.pressure,
+        )
+        if isinstance(mass_rate, (int, float)):
+            mass_rate_array = np.array([mass_rate])
+        else:
+            mass_rate_array = mass_rate
+        mass_rate_array = cond._apply_skip_partitioning(mass_rate_array)
+        norm_conc = particle_data.concentration[0] / particle_data.volume[0]
+        mass_transfer = get_mass_transfer(
+            mass_rate=mass_rate_array,
+            time_step=1.0,
+            gas_mass=np.asarray(gas_data.concentration[0], dtype=np.float64),
+            particle_mass=particle_data.masses[0],
+            particle_concentration=norm_conc,
+        )
+        species_mass = particle_data.masses[0]
+        species_count = 1 if species_mass.ndim == 1 else species_mass.shape[1]
+        if mass_transfer.ndim == 1:
+            mass_transfer = mass_transfer.reshape(-1, species_count)
+        elif mass_transfer.shape[1] > species_count:
+            mass_transfer = mass_transfer[:, :species_count]
+        elif mass_transfer.shape[1] < species_count:
+            mass_transfer = np.broadcast_to(
+                mass_transfer, (mass_transfer.shape[0], species_count)
+            )
+
+        cond.step(
+            particle=particle_data,
+            gas_species=gas_data,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+
+        latent_heat = cond._latent_heat_strategy.latent_heat(self.temperature)
+        np.testing.assert_allclose(
+            cond.last_latent_heat_energy,
+            np.sum(mass_transfer * latent_heat),
+            rtol=1e-14,
+        )
+
+    def test_step_energy_tracking_skip_partitioning(self):
+        """Skip-partitioning affects latent heat bookkeeping."""
+        particle = copy.deepcopy(self.particle)
+        gas_species = copy.deepcopy(self.gas_species)
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+            skip_partitioning_indices=[1],
+        )
+
+        particle_data, _ = _unwrap_particle(particle)
+        gas_data, _ = _unwrap_gas(gas_species)
+        mass_rate = cond.mass_transfer_rate(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+        )
+        if isinstance(mass_rate, (int, float)):
+            mass_rate_array = np.array([mass_rate])
+        else:
+            mass_rate_array = mass_rate
+        mass_rate_array = cond._apply_skip_partitioning(mass_rate_array)
+        norm_conc = particle_data.concentration[0] / particle_data.volume[0]
+        mass_transfer = get_mass_transfer(
+            mass_rate=mass_rate_array,
+            time_step=1.0,
+            gas_mass=np.asarray(gas_data.concentration[0], dtype=np.float64),
+            particle_mass=particle_data.masses[0],
+            particle_concentration=norm_conc,
+        )
+        species_mass = particle_data.masses[0]
+        species_count = 1 if species_mass.ndim == 1 else species_mass.shape[1]
+        if mass_transfer.ndim == 1:
+            mass_transfer = mass_transfer.reshape(-1, species_count)
+        elif mass_transfer.shape[1] > species_count:
+            mass_transfer = mass_transfer[:, :species_count]
+        elif mass_transfer.shape[1] < species_count:
+            mass_transfer = np.broadcast_to(
+                mass_transfer, (mass_transfer.shape[0], species_count)
+            )
+
+        cond.step(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+
+        latent_heat = cond._latent_heat_strategy.latent_heat(self.temperature)
+        np.testing.assert_allclose(
+            cond.last_latent_heat_energy,
+            np.sum(mass_transfer * latent_heat),
+            rtol=1e-14,
+        )
+
+    def test_step_energy_sign_convention(self):
+        """Condensation is positive energy; evaporation is negative."""
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+        )
+        gas_species, particle, _, _ = self._build_single_species_fixture()
+
+        cond.step(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+        self.assertGreater(cond.last_latent_heat_energy, 0.0)
+
+        vp_water = par.gas.VaporPressureFactory().get_strategy("water_buck")
+        mm_water = 18.015e-3
+        gas_sub = (
+            par.gas.GasSpeciesBuilder()
+            .set_name("H2O")
+            .set_molar_mass(np.array([mm_water]), "kg/mol")
+            .set_vapor_pressure_strategy([vp_water])
+            .set_concentration(np.array([1e-30]), "kg/m^3")
+            .set_partitioning(True)
+            .build()
+        )
+        particle_evap = copy.deepcopy(particle)
+        cond.step(
+            particle=particle_evap,
+            gas_species=gas_sub,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+        self.assertLess(cond.last_latent_heat_energy, 0.0)
+
+    def test_step_isothermal_parity(self):
+        """Latent heat fallback matches isothermal step outputs."""
+        particle_latent = copy.deepcopy(self.particle)
+        gas_latent = copy.deepcopy(self.gas_species)
+        particle_iso = copy.deepcopy(self.particle)
+        gas_iso = copy.deepcopy(self.gas_species)
+        latent_cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+        )
+        iso_cond = CondensationIsothermal(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+        )
+
+        particle_out, gas_out = latent_cond.step(
+            particle=particle_latent,
+            gas_species=gas_latent,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+        particle_iso_out, gas_iso_out = iso_cond.step(
+            particle=particle_iso,
+            gas_species=gas_iso,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+
+        np.testing.assert_allclose(
+            particle_out.get_mass(), particle_iso_out.get_mass(), rtol=1e-15
+        )
+        np.testing.assert_allclose(
+            gas_out.get_concentration(),
+            gas_iso_out.get_concentration(),
+            rtol=1e-15,
+        )
+        self.assertEqual(latent_cond.last_latent_heat_energy, 0.0)
+
+    def test_step_returns_correct_types_legacy(self):
+        """Legacy inputs return legacy types."""
+        cond = CondensationLatentHeat(molar_mass=self.molar_mass)
+
+        particle_new, gas_new = cond.step(
+            particle=copy.deepcopy(self.particle),
+            gas_species=copy.deepcopy(self.gas_species),
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+
+        self.assertIsInstance(particle_new, ParticleRepresentation)
+        self.assertIsInstance(gas_new, GasSpecies)
+
+    def test_step_returns_correct_types_data(self):
+        """Data-only inputs return data types."""
+        particle_data, gas_data = self._make_data_inputs()
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            activity_strategy=self.activity_strategy,
+            surface_strategy=self.surface_strategy,
+            vapor_pressure_strategy=self.vapor_pressure_strategy,
+        )
+
+        particle_new, gas_new = cond.step(
+            particle=particle_data,
+            gas_species=gas_data,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+
+        self.assertIsInstance(particle_new, ParticleData)
+        self.assertIsInstance(gas_new, GasData)
+
+    def test_step_updates_gas_when_enabled(self):
+        """Gas concentrations decrease by mass transfer when enabled."""
+        particle = copy.deepcopy(self.particle)
+        gas_species = copy.deepcopy(self.gas_species)
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+            update_gases=True,
+        )
+
+        particle_data, _ = _unwrap_particle(particle)
+        gas_data, _ = _unwrap_gas(gas_species)
+        mass_rate = cond.mass_transfer_rate(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+        )
+        mass_rate_array = (
+            np.array([mass_rate])
+            if isinstance(mass_rate, (int, float))
+            else mass_rate
+        )
+        mass_rate_array = cond._apply_skip_partitioning(mass_rate_array)
+        norm_conc = particle_data.concentration[0] / particle_data.volume[0]
+        expected_mass_transfer = get_mass_transfer(
+            mass_rate=mass_rate_array,
+            time_step=1.0,
+            gas_mass=np.asarray(gas_data.concentration[0], dtype=np.float64),
+            particle_mass=particle_data.masses[0],
+            particle_concentration=norm_conc,
+        )
+        expected_gas = (
+            gas_species.get_concentration() - expected_mass_transfer.sum(axis=0)
+        )
+
+        _, gas_new = cond.step(
+            particle=particle,
+            gas_species=gas_species,
+            temperature=self.temperature,
+            pressure=self.pressure,
+            time_step=1.0,
+        )
+
+        np.testing.assert_allclose(
+            gas_new.get_concentration(), expected_gas, rtol=1e-15
+        )
+
+    def test_step_clamps_gas_concentration_nonnegative(self):
+        """Gas concentrations are clamped to nonnegative values."""
+        particle_data, gas_data = self._make_data_inputs()
+        gas_data.concentration[0] = np.array([1e-12, 2e-12])
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            activity_strategy=self.activity_strategy,
+            surface_strategy=self.surface_strategy,
+            vapor_pressure_strategy=self.vapor_pressure_strategy,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+        )
+        mass_transfer = np.full_like(particle_data.masses[0], 1.0)
+
+        with patch(
+            "particula.dynamics.condensation.condensation_strategies."
+            "get_mass_transfer",
+            return_value=mass_transfer,
+        ):
+            _, gas_new = cond.step(
+                particle=particle_data,
+                gas_species=gas_data,
+                temperature=self.temperature,
+                pressure=self.pressure,
+                time_step=1.0,
+            )
+
+        np.testing.assert_allclose(gas_new.concentration[0], 0.0, rtol=1e-15)
+
+    def test_step_dynamic_viscosity_pass_through(self):
+        """Dynamic viscosity override is forwarded to mass_transfer_rate."""
+        particle = copy.deepcopy(self.particle)
+        gas_species = copy.deepcopy(self.gas_species)
+        cond = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+        )
+        override = 1.23e-5
+
+        with patch.object(
+            cond, "mass_transfer_rate", wraps=cond.mass_transfer_rate
+        ) as spy:
+            cond.step(
+                particle=particle,
+                gas_species=gas_species,
+                temperature=self.temperature,
+                pressure=self.pressure,
+                time_step=1.0,
+                dynamic_viscosity=override,
+            )
+
+        self.assertEqual(spy.call_args.kwargs["dynamic_viscosity"], override)
+
+    def test_step_mass_condensation_runnable_compatibility(self):
+        """MassCondensation executes with latent heat strategies."""
+        particle = copy.deepcopy(self.particle)
+        gas_species = copy.deepcopy(self.gas_species)
+        atmosphere = par.gas.Atmosphere(
+            temperature=self.temperature,
+            total_pressure=self.pressure,
+            partitioning_species=gas_species,
+            gas_only_species=[],
+        )
+        aerosol = par.Aerosol(atmosphere=atmosphere, particles=particle)
+
+        cond_latent = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+            latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6),
+        )
+        condensation = par.dynamics.MassCondensation(
+            condensation_strategy=cond_latent
+        )
+        aerosol_out = condensation.execute(aerosol, time_step=1.0)
+
+        self.assertIsInstance(aerosol_out.particles, ParticleRepresentation)
+        self.assertIsInstance(
+            aerosol_out.atmosphere.partitioning_species, GasSpecies
+        )
+
+        aerosol_iso = par.Aerosol(
+            atmosphere=copy.deepcopy(atmosphere),
+            particles=copy.deepcopy(self.particle),
+        )
+        cond_iso = CondensationLatentHeat(
+            molar_mass=self.molar_mass,
+            diffusion_coefficient=self.diffusion_coefficient,
+            accommodation_coefficient=self.accommodation_coefficient,
+        )
+        condensation_iso = par.dynamics.MassCondensation(
+            condensation_strategy=cond_iso
+        )
+        aerosol_iso_out = condensation_iso.execute(aerosol_iso, time_step=1.0)
+
+        self.assertIsInstance(aerosol_iso_out.particles, ParticleRepresentation)
+        self.assertIsInstance(
+            aerosol_iso_out.atmosphere.partitioning_species, GasSpecies
+        )
diff --git a/readme.md b/readme.md
index dacb07973..e08978f18 100644
--- a/readme.md
+++ b/readme.md
@@ -99,7 +99,8 @@ particula/
   `particula.dynamics.get_latent_heat_energy_released`
 - **Condensation Strategies** — `CondensationIsothermal` plus
   `CondensationLatentHeat` with latent-heat-corrected
-  `mass_transfer_rate()`/`rate()` (step integration pending)
+  `mass_transfer_rate()`/`rate()` and `step()` energy tracking via
+  `last_latent_heat_energy`, with optional `dynamic_viscosity` override
 - **Latent Heat Factories** — Build constant, linear, and power-law latent heat
   strategies via `particula.gas.LatentHeatFactory` with unit-aware builders and
   gas-phase exports for upcoming non-isothermal workflows